ORIGINAL RESEARCH

published: 14 February 2019

doi: 10.3389/frobt.2019.00007

Attitude Stabilization of Spacecraft

in Very Low Earth Orbit by
Center-Of-Mass Shifting
Josep Virgili-Llop, Halis C. Polat and Marcello Romano*

Spacecraft Robotics Laboratory, Naval Postgraduate School, Monterey, CA, United States

At very low orbital altitudes (.450 km) the aerodynamic forces can become major
attitude disturbances. Certain missions that would benefit from a very low operational
altitude require stable attitudes. The use of internal shifting masses, actively shifting the
location of the spacecraft center-of-mass, thus modulating, in direction and magnitude,
the aerodynamic torques, is here proposed as a method to reject these aerodynamic
disturbances. A reduced one degree-of-freedom model is first used to evaluate the
disturbance rejection capabilities of the method with respect to multiple system
parameters (shifting mass, shifting range, vehicle size, and altitude). This analysis shows
that small shifting masses and limited shifting ranges suffice if the nominal center-of-mass
is relatively close to the estimated center-of-pressure. These results are confirmed when
the analysis is extended to a full three rotational degrees-of-freedom model. The use
of a quaternion feedback controller to detumble a spacecraft operating at very low
altitudes is also explored. The analysis and numerical simulations are conducted using a
Edited by:
Riccardo Bevilacqua,
nonlinear dynamic model that includes the full effects of the shifting masses, a realistic
University of Florida, United States atmospheric model, and uncertain spacecraft aerodynamic properties. Finally, a practical
Reviewed by: implementation on a 3U CubeSat using commercial-off-the-shelf components is briefly
Marco Sabatini,
presented, demonstrating the implementation feasibility of the proposed method.
Sapienza University of Rome, Italy
David Andres Perez, Keywords: spacecraft aerodynamics, attitude stabilization, Very Low Earth Orbit, attitude control, shifting masses,
General Motors, United States movable masses, CubeSat, aerodynamic disturbance
*Correspondence:
Marcello Romano
mromano@nps.edu 1. INTRODUCTION
Specialty section: Lowering the operational altitude of Earth observation spacecraft can increase the overall
This article was submitted to cost-effectiveness of a space system (Shao et al., 2014). For example, by lowering the operational
Space Robotics, altitude, the resolution of a given optical instrument, its radiometric performance and the
a section of the journal geospatial accuracy of its imagery are improved. For radar payloads, either the antenna size or
Frontiers in Robotics and AI
the transmission power can be reduced. Furthermore, a given launcher can usually deliver more
Received: 14 June 2018 payload at lower altitude orbits or, for a given spacecraft mass, a less capable and potentially more
Accepted: 21 January 2019 cost-effective launcher can be used (Virgili-Llop, 2014; Virgili-Llop et al., 2014a).
Published: 14 February 2019
Lowering the operational altitude forces spacecraft to orbit through denser regions of the
Citation: atmosphere. The interaction of spacecraft with the residual atmosphere results in aerodynamic
Virgili-Llop J, Polat HC and
forces that, at low altitudes, can become major orbit and attitude perturbations (Fortescue and
Romano M (2019) Attitude
Stabilization of Spacecraft in Very Low
Stark, 1995). As the aerodynamic forces are only dominant in the lower part of the Low Earth Orbit
Earth Orbit by Center-Of-Mass (LEO) range the term Very Low Earth Orbit (VLEO) is used in this paper to make clear that the
Shifting. Front. Robot. AI 6:7. considered orbit range extends only up to ∼450 km in altitude (Virgili-Llop, 2014; Virgili-Llop
doi: 10.3389/frobt.2019.00007 et al., 2014a).

Frontiers in Robotics and AI | www.frontiersin.org 1 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

Although these aerodynamic perturbations are usually Of particular interest is the work by Chesi et al. (2017) who
perceived as drawbacks, they can also be seen as an advantage proposes the use of aerodynamic drag to generate attitude control
(Virgili-Llop, 2014; Virgili-Llop et al., 2014a). As the orbital torques modulated in magnitude and direction by actively
lifetime is reduced, there is no need to de-orbit spacecraft at shifting a set of internal masses. Although Chesi’s work, simplifies
the end of their operational life. Space debris also decay at a the effects of the shifting masses on the spacecraft dynamics,
faster rate and thus do not accumulate at the same pace in ignores the variable and unpredictable nature of the Earth’s
VLEO, reducing the collision risk and greatly increasing the atmosphere, and assumes that the aerodynamic properties are
required object density to generate a Kessler syndrome runaway known and constant, it shows the conceptual feasibility of using
(Wertz et al., 2012). Additionally, the aerodynamic forces can shifting masses to control the aerodynamic torques. In particular,
be harnessed to control the spacecraft’s orbit (Bevilacqua and it shows that by using a set of three shifting masses augmented
Romano, 2008; Virgili-Llop et al., 2014b; Virgili et al., 2015) and by reaction wheels or magnetic torquers and using an adaptive
attitude (Kumar et al., 1995, 1996; Psiaki, 2004; Gargasz, 2007; non-linear feedback control law, a spacecraft could be slowly
Guettler, 2007). brought, from any initial attitude and angular velocity, to a
Many missions that would benefit from lower operational desired attitude whilst minimizing the use of the reaction wheels
altitudes also require a constant and stable attitude (e.g., geodesy or magnetic torquers.
spacecraft, Earth observation). In such missions, the attitude The work presented in this paper takes Chesi’s concept one
perturbations caused by the aerodynamic forces need to be step forward by dropping the dynamic model simplifications,
eliminated. An effective and conceptually simple measure to introducing uncertainties into the aerodynamic properties, and
reduce the aerodynamic forces for a given operational altitude adding atmospheric variability. Additional contributions of the
is to minimize the spacecraft’s cross section area exposed to the work presented in this paper are a sensitivity analysis of the
incident flow. For this very reason, spacecraft operating in VLEO method’s performance with respect to the CoP to CoM distance,
tend to be slender (Bowman and Lewis, 2002; Drinkwater et al., size of the spacecraft, and operating altitude. Additionally, an
2007). To minimize the attitude perturbation it is also highly assessment of the implementation feasibility of the concept is
desirable to design spacecraft with their Center-of-Mass (CoM) briefly presented. As in Chesi et al. (2017), the assumptions used
as close as possible to the spacecraft’s Center-of-Pressure (CoP), to derive the dynamic model and the controllers are made explicit
thus minimizing the force lever arm and ultimately reducing throughout the paper and are marked with Asm.
the aerodynamic disturbance torque. Unfortunately, spacecraft This paper is organized as follows. The spacecraft model is
aerodynamics uncertainties (Moe and Moe, 2010; Prieto et al., briefly presented in section 2. Then the equations of motion of
2014) and atmospheric variability (Larsen and Fesen, 2009; a spacecraft with internal moving parts are derived in section
Pardini et al., 2012) introduce uncertainties and variability to 3. The uncertain nature of the aerodynamic disturbance caused
the CoP location. Additional practical design constrains on the by a variable atmosphere and the uncertain aerodynamics is
location of the CoM make the realization of an overlapping subsequently presented in section 4. A reduced one rotational
CoP and CoM impossible in practice. A residual aerodynamic degree-of-freedom model with one shifting mass driven by a
disturbance torque will remain and will need to be rejected. Proportional-Integral-Derivative (PID) controller is derived in
These attitude perturbations can be compensated for by using section 5. This PID controller is used to analyze the disturbance
traditional attitude control actuators (such as reaction wheels). rejection capabilities of the system with respect to several
At VLEO the aerodynamic disturbances can be significant and parameters (shifting mass, shifting range, operating altitude and
can present a secular component that can quickly saturate vehicle size). Then we use a full three degrees-of-freedom model
momentum exchange devices. In this paper, we explore the use with two shifting masses driven by a Linear Quadratic Regulator
of internal shifting masses as a method to control and ultimately (LQR) based controller moving along the pitch and yaw axes
reject these undesired aerodynamic disturbances. The set of and augmented by an ideal actuator in roll in section 6 to
internal shifting internal masses, actively change the spacecraft’s confirm that the results obtained in the one rotational degree-of-
CoM location, modulating, in direction and magnitude, the freedom reduced model also apply in a three degrees-of-freedom
aerodynamic torques. Specifically, we are interested in: devising model. Then, the detumbling capabilities of the proposed
control methods to drive the shifting masses, evaluating the method are briefly explored in section 6 with a quaternion
disturbance rejection capabilities under realistic conditions, and feedback controller. Finally, a practical implementation, only
evaluating the implementation feasibility of the whole shifting using Commercial-Off-The-Shelf (COTS) components, on a 3U
masses concept. CubeSat is presented in section 7.
The use of shifting masses as attitude control actuators has
already been proposed in the past to help detumble spacecraft
(Edwards and Kaplan, 1974; Kunciw and Kaplan, 1976), control 2. SPACECRAFT MODEL
the coning motion of a spinning spacecraft (Hamidi-Hashemi,
1993; Halsmer and Mingori, 1995; Janssens and van der Ha, To keep the analysis as general as possible, a spherically shaped
2014), control the pitch and yaw of solar-sails (Wie, 2004; Wie spacecraft has been used. Although a spheric spacecraft may be
and Murphy, 2007; Scholz et al., 2011) and, in general, to initially perceived as a simplistic case, it can already be used
complement traditional attitude control actuators (Kumar, 2010; to illustrate the effects of aerodynamic uncertainties without
Ahn, 2012; Atkins and Henderson, 2012). dwelling into more complex shapes. Also, the simple relationship

Frontiers in Robotics and AI | www.frontiersin.org 2 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

between the size of the sphere (its radius) and its mass properties atmospheric co-rotation). The implications of this assumption
(inertia and mass) is used to extract the general trends with will become clear in section 4.4, but briefly stated, by selecting an
respect to spacecraft size. aerodynamic equilibrium attitude we avoid secular aerodynamic
torques. If an arbitrary attitude was selected the location of the
• Asm.1: The spacecraft is spherically shaped.
shifting masses would be biased in order to provide this secular
The spherically shaped spacecraft hosting the shifting masses torque (essentially moving the system’s CoM and forcing the
(host spacecraft) is assumed to be composed of a homogeneous desired attitude to become an equilibrium one).
density sphere and a fixed discrete point mass (not a shifting
• Asm.2: The desired attitude is an aerodynamic equilibrium
mass) as shown in Figure 1. The discrete point mass is added
attitude (in the absence of wind and atmospheric co-rotation).
to the host vehicle to obtain a host spacecraft CoM that is
not coincident with the sphere’s geometric center. It is worth The inertia properties of the host spacecraft J 0 can be computed
mentioning at this point that the shifting masses are not part using Equation (1), with the distance between the point mass and
of the host spacecraft. Excluding the shifting masses greatly the sphere’s center denoted by dMP . If dMP > 0 the CoM will be
simplifies the equations of motion as it is shown in section 2. located in the positive side of î and if dMP < 0 then the CoM will
The mass of a homogeneous density sphere is MS = be in the negative side of î. The location of the host spacecraft’s
ρS 4/3πR3 and its inertia J S = Iρs 8/15πR5 , with ρS denoting the CoM from the sphere center is simply dCoM0 =
κdMP
(1+κ) .
sphere’s density, R the sphere’s radius, and I the identity matrix.
The mass of the fixed discrete point mass MP can be expressed
0 0 0
 
as a fraction κ of the homogeneous sphere’s mass MP = κMS . MS (1 + κ)  2
J0 = JS + 0 dMP 0  (1)
The total mass of the host vehicle M0 is simply M0 = MS (1 + κ). κ 2
0 0 dM
Let the host spacecraft’s CoM define the origin of the spacecraft’s p

body frame B0 with the discrete point mass located along the î
axis, as depicted in Figure 1. Note that the host spacecraft is symmetric with respect to the
The goal of the attitude control system is to keep the attitude roll axis î and thus the definition of the pitch ĵ and yaw k̂ axes
of the spacecraft stable with respect to the orbital frame (k̂orbit is arbitrary.
points nadir, îorbit along the inertial velocity vector, and ĵorbit To this host vehicle, whose mass and inertia properties are
completes the right hand triad). In this case, the desired attitude fixed, known and constant, a set of N shifting masses mn is added,
will be to align the axes of the body frame B0 with the axes of the altering the mass, inertia, and CoM location of the resulting
orbital frame, thus in the desired attitude î represents the roll axis, combined system. As the body reference frame is defined with
respect to the CoM of the host spacecraft, excluding the shifting
ĵ the pitch axis, and k̂ the yaw axis.
masses, the system’s CoM, including the shifting masses, will be
As the location of the discrete point mass is artificially
variable and thus will not be located at the origin of the body
restricted to the î axis, the desired attitude represents an
reference frame B0 .
aerodynamic equilibrium attitude (in the absence of wind and

3. DYNAMIC MODEL
The equations of motion for a system of connected rigid bodies
is derived in this section. These equations of motion, which
take into account the dynamic effects of the shifting masses,
are used during the numerical simulations and serve as the
starting dynamic model used to derive the controllers. The
following assumptions are made to simplify the formulation of
the dynamics.
• Asm.3: The host spacecraft is a rigid body.
• Asm.4: The shifting masses are rigid bodies.
Under these assumptions, the fundamental equation describing
the rotational motion of a system of connected rigid bodies is
given by Equation (2) (Grubin, 1962; Edwards and Kaplan, 1974).

τ = Ḣ + S × a (2)

In Equation (2) τ denotes the external torques around an

arbitrary reference point, Ḣ the time derivative of the system’s
angular momentum around this reference point, S the system’s
FIGURE 1 | Spacecraft model with its homogeneous sphere, the discrete
point mass along the î axis and a generic shifting mass mn .
first moment of mass with respect to the reference point, and
a denotes the inertial acceleration of the reference point. The

Frontiers in Robotics and AI | www.frontiersin.org 3 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

system’s reference point can be arbitrarily selected and can move and shifting masses angular momentum time derivatives (using
in an arbitrary manner. It is interesting to note that if the the transport theorem where appropriate).
acceleration of the reference point is a = 0 or if the reference
point is selected as the system’s CoM (S = 0) then the usual N
X
expression τ = Ḣ is recovered. Ḣ = Ḣ 0 + Ḣ n (8)
As introduced in the previous section, the system is composed n=1
of a host vehicle to which N shifting masses are added. To
simplify the derivation of the equations, and without loosing Ḣ 0 = J 0 ω̇0 + ω0 × H 0 (9)
generality, it will be considered that the reference point of the
system is the host vehicle’s CoM (which excludes the shifting
masses). This reference point is the origin of the body frame B0 Ḣ n = J n ω̇n + ωn × H n + mn rn × r̈n (10)
as shown in Figure 1. This assumption is useful, as by definition,
The inertial acceleration of the shifting masses can be computed
the host vehicle’s mass M0 and its inertia J 0 are constant when
with the transport theorem, resulting in the following well
projected into B0 . Additionally, the movement of the shifting
known equation.
masses can be easily known with respect to the host spacecraft
reference frame B0 .
r̈n = ω0 × (ω0 × rn ) + ω̇0 × rn + 2ω0 × ṙ′n + r̈′n (11)
The N shifting masses have their own reference frames Bn with
their origin located at the CoM of the shifting mass. The shifting
Note how ṙ′n and r̈′n are the shifting masses relative velocity and
masses can be rigid bodies or point masses. If they are rigid bodies
acceleration with respect to B0 . These ṙ′n and r̈′n magnitudes can
the inertia of the shifting masses in Bn will be denoted as J Bn n and be measured inside the host spacecraft.
when expressed in B0 it will be simply denoted as J n . Moving on with the other terms in Equation (2), the first
The location of a shifting mass with respect to the reference moment of mass is defined as in the following equation.
point (the origin of B0 ) will be denoted as rn . The ṙn and r̈n
terms will denote the inertial linear velocity and acceleration of N
X
the shifting mass expressed in the B0 frame. The inertial angular S= mn r n (12)
velocity of the host vehicle frame B0 is denoted by ω0 and the n=1
term ωn denotes the inertial angular velocities of the Bn frames
expressed in the B0 frame. The inertial acceleration of the origin of B0 (the reference point
The inertial angular velocity of the shifting mass can be in Equation 2) can then be written as follows.
computed using the following equation, with ω′n being the
relative angular velocity of the shifting mass reference Bn with a = r̈B0 = r̈CoM − r̈′CoM (13)
respect to the host vehicle reference B0 .
The r̈′CoM term denotes the acceleration of the system’s CoM with
ωn = ω0 + ω′n (3) respect to B0 (the relative movement of the system’s CoM) and
r̈CoM is the inertial acceleration of the system’s CoM (due to the
In such a multibody system, the total angular momentum H is external forces F). The r̈CoM acceleration can be easily computed
composed of the sum of the angular momentum of host vehicle using Newton’s second law and r̈′CoM is obtained by computing
H 0 and of the shifting masses H n , as shown in Equation (4). the relative CoM acceleration as follows.

N F
X r̈CoM = PN (14)
H = H0 + Hn (4) M0 + n=1 mn
n=1

PN
n=1 mn r̈ n
H 0 = J 0 ω0 (5) r̈′CoM = (15)
M0 + N
P
n=1 mn

H n = J n ωn + mn rn × ṙn (6) With the equations above, Equation (2) can be fully expanded as
in Equation (16).
The linear inertial velocity of the shifting mass ṙn can be simply
computed using the transport theorem, resulting in the following N N
equation, where ṙ′n denotes the relative velocity of the shifting
X X
J 0 ω̇0 + ω0 × J 0 ω0 + J n ω̇n + ωn × (J n ωn + mn rn × ṙn )
mass with respect to B0 . n=1 n=1
N N
!
ṙn = ṙ′n + ω0 × rn (7)
X 1 X
+ (mn rn × r̈n ) + ...... + PN mn r̈n
n=1 M0 + n=1 mn n=1
To use Equation (2) the angular momentum needs to be N N
differentiated. Deriving (Equation 4) it follows that the total
X F X
× mn r n = τ + PN × mn r n (16)
angular momentum time derivative is the sum of the host vehicle n=1 M0 + n=1 mn n=1

Frontiers in Robotics and AI | www.frontiersin.org 4 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

It is interesting to note that (Equation 16) also contains the generic force coefficients Cf , the drag CD (anti velocity), and
case of momentum exchange devices (reaction wheels and lift CL (normal to velocity) coefficients, which leads to the well
control moment gyroscopes) and thus these devices can be easily known drag and lift (Equations 21, 22).
incorporated into this analysis.
The aerodynamic effects on the attitude dynamics of the 1 2
F= ρV Aref Cf (20)
system are represented by the external torques τ and forces F. 2 ∞
It is important to note that the torques τ are computed around
the fixed reference point (the host vehicle CoM) and not with 1 2
respect to the moving system’s CoM. The term in Equation (16) D= ρV Aref CD (21)
2 ∞
that contains the external forces F accommodates this difference.
It is also important to note that (Equation 16) contains several
1 2
terms that depend on the shifting masses relative velocities and L= ρV Aref CL (22)
accelerations, thus accommodating the dynamic effects of the 2 ∞
shifting masses. From Equation 20 (or Equations 21, 22) the atmospheric density
ρ, the relative velocity with respect to the flow V∞ and the
3.1. Point Mass Simplification force coefficients CD and CL need to be estimated (using the
A very useful simplification is obtained when it is assumed that environment and gas-surface interaction models) before the
the shifting masses are point masses and do not posses inertia aerodynamic forces can be computed.
J n = 0. In such a case, the rotation of the shifting masses is It is worth noting at this point that the atmospheric density
irrelevant and their translation is the only parameter that affect approximately increases exponentially with decreasing altitude
the dynamics. Under such assumption the general equations of and thus the aerodynamic forces magnitude will also increase
motion (Equation 16) can be simplified as in Equation (17). exponentially with decreasing orbital altitude. The aerodynamic
disturbance is therefore strongly dependent on the altitude and
N
X dominates at very low orbital altitudes.
J 0 ω̇0 + ω0 × J 0 ω0 + (mn rn × r̈n )
Ideally, the controller that regulates the shifting masses
n=1
position would know the direction and magnitude of the relative
N N
!
1 X X flow V∞ , the atmospheric density ρ, and the aerodynamic
+ mn r̈n × mn rn = ......
M0 +
PN properties of the spacecraft Cf . With this information it would
n=1 mn n=1 n=1
be able to accurately estimate the aerodynamic torque that
N
F X the spacecraft is experiencing and drive the shifting masses
= τ+ × mn r n (17)
M0 +
PN to reject it. Unfortunately, the atmospheric environment is
n=1 mn n=1
highly variable and poorly predictable (Larsen and Fesen, 2009;
For a single point mass and introducing the concept of reduced Pardini et al., 2012) and spacecraft aerodynamics are not
mass µ (Equation 18), the equation can be further simplified particularly well understood (Moe and Moe, 2010; Prieto et al.,
to finally obtain (Equation 19), recovering the expression from 2014). As a consequence, the controller will not be able to
Edwards and Kaplan (1974). obtain accurate estimates of the aerodynamic torque magnitude
or direction.
mM0 In the numerical simulations conducted in this paper, state
µ= (18) of the art atmospheric and spacecraft aerodynamic models (see
M0 + m
following sections) are used to obtain what it is assumed to be the
truth values. The controller will then estimate these magnitudes
µF using simplified aerodynamics and atmospheric models. This set-
J 0 ω̇0 + ω0 × J 0 ω0 + µr × r̈ = τ + ×r (19)
M0 up ensures the presence of realistic atmospheric variability and
realistic aerodynamic properties, while emulating the uncertainty
4. AERODYNAMIC MODELING that a controller will be subjected to.

The residual atmosphere present at orbital altitudes causes 4.1. Atmospheric Density Model
spacecraft to experience aerodynamic forces (mainly For this study, the NRLMSISE-00 (Picone et al., 2002)
aerodynamic drag). Orbital decay is the main effect of atmospheric model is used to estimate the atmospheric density
aerodynamic drag but these aerodynamic forces will also induce ρ. This specific atmospheric model offers a good balance
aerodynamic torques and thus perturb the spacecraft’s attitude. between model accuracy and computational complexity (ECSS
In general, Equation (20) can be used to compute these Secretariat, 2008).
aerodynamic forces, where ρ denotes the atmospheric density,
• Asm.5: The atmosphere density behaves as predicted by the
V∞ the relative velocity of the spacecraft with respect to the flow,
NRLMSISE-00 model.
Aref an arbitrary reference area (usually taken as the spacecraft’s
cross section area), and Cf the force coefficients (along the three The Earth’s atmosphere not only exhibits vertical density
different axes). Special cases of Equation (20) use, instead of the variations but also horizontal ones (as the day-to-night density

Frontiers in Robotics and AI | www.frontiersin.org 5 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

changes among others). Thus, a spacecraft orbiting in a circular • Asm.7: The atmosphere co-rotates and its wind behaves as
orbit will experience density variations (that will modify the predicted by the HWM07 model.
magnitude of the aerodynamic forces). Figure 2 shows the
Figure 3B shows the sideslip angle caused by the atmospheric
density variations, with respect to the orbit’s mean density
co-rotation and wind (using the HWM07 model) assuming that
(using the NRLMSISE-00 model), for a 10:30 Local Time
a spacecraft is aligned with its inertial velocity in a 10:30 LTAN
of the Ascending Node (LTAN) circular Sun-synchronous
circular Sun-synchronous orbit at different altitudes in moderate
orbit at different orbital altitudes in moderate solar activity
solar activity (ISO 14222, 2013).
(ISO 14222, 2013). Figure 2 exemplifies how variable the density
is and thus how variable the magnitude of the drag, and
consequently the aerodynamic torque, is during these typical
Sun-synchronous orbits. 4.3. Gas-Surface Interaction Model
In the orbital environment (>200 km in altitude) the residual
4.2. Wind Model atmosphere can no longer be considered as a continuum but,
It is not uncommon to assume that a spacecraft’s inertial velocity given its low density, needs to be considered as a rarefied-gas
is equal to the spacecraft’s relative velocity with respect to the (Bird, 1994). The mean free path λ of an atmospheric gas particle
incoming flow. This assumption ignores that the atmosphere co- is, in general, much greater than a representative spacecraft
rotates with the Earth (Challinor, 1968; King-Hele, 1987, 1992) dimension (λ > 100 m at 200 km altitude; Virgili-Llop, 2014).
and that there is atmospheric time-varying wind (Killeen et al., Consequently, it can be assumed that the interactions between
1982; King-Hele and Walker, 1988). These two effects will make gas particles (collisions) are very rare, and thus they can be
the direction and magnitude of the relative flow V∞ differ, in safely neglected. Therefore, the Gas-Surface Interactions (GSI)
direction and magnitude, from the inertial velocity. completely dominate the interaction of the spacecraft with its
The atmospheric wind is also highly variable, spatially and surrounding gas.
temporally. Figure 3A shows an example distribution of the
• Asm.8: Gas-gas particle interactions are negligible.
wind. As the atmospheric wind has not been as extensively
studied as other atmospheric properties, the existing models The GSI are dependent on several gas and surface parameters.
are less accurate (Larsen and Fesen, 2009). In this work, the As these interactions occur at the molecular scale, molecular
HWM07 (Drob et al., 2008) wind model is used. It has to scale parameters are also relevant (e.g., lattice configuration and
be noted that this model only provides zonal and meridional surface roughness among others). The high thermal velocity
wind profiles representative of the climatological averages for of the gas particles (∼1 km/s at 350 km), due to the high
various geophysical conditions. Vertical winds, which usually temperature and low density of the gas, produces a flow that is not
have smaller magnitudes, are not included in the model. Real collimated. The non collimated flow leads to particles colliding
instantaneous values may show finer temporal and spatial with surfaces that would intuitively appear to be shadowed from
variations than the ones provided by the model and their effects the flow.
would need to be considered if this concept is brought to There are several GSI models (Bird, 1994) and in this study
operational use. the Sentman model (Sentman, 1961) will be used as it is the de
facto standard to compute spacecraft aerodynamic coefficients
at low altitudes (Moe and Moe, 2005, 2010). A comprehensive
description of the models used in spacecraft aerodynamics can
be found elsewhere (Moe and Moe, 2010; Prieto et al., 2014).
• Asm.9: The lift and drag coefficients behave according to the
Sentman model.
The Sentman model takes into account the thermal velocity
distribution of the gas particles and assumes that all the incident
gas particles that collide with a surface are adsorbed to be later
diffusely reemitted. In the LEO range this seems to be true from
the limited available orbital data (Gregory and Peters, 1987; Moe
et al., 1998). The particles are then reemitted with partial thermal
equilibrium with the spacecraft surface. The degree of thermal
equilibrium is denoted by the energy accommodation coefficient
σa . In this model, the Cd and Cl can be written, following a
notation similar to Sutton (2009) and Doornbos (2011), as in
Equations (23, 25).

FIGURE 2 | Density variations during a typical circular Sun-synchronous orbit

√
 
at different operational altitudes.
P γ vre

Cd = √ + γ QZ + γ πZ + P (23)
π 2 V∞

Frontiers in Robotics and AI | www.frontiersin.org 6 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

FIGURE 3 | Wind pattern according to the HWM93 (Hedin et al., 1996) model at 450 km with moderate solar activity (ISO 14222, 2013) during northern hemisphere
summer solstice (A) and yaw angle caused by atmospheric co-rotation and wind during a circular sun-synchronous orbit at different operational altitudes (B).

R
Cd dA In the VLEO range the atomic oxygen is one of the dominant
CD = (24)
Aref species. These atomic oxygen gas molecules get adsorbed into
the spacecraft surfaces masking the original surface properties.
√
 
l vre
Having a surface covered with atomic oxygen rises the
Cl = lGZ + γ πZ + P (25)
2 V∞ accommodation coefficient to a level between 0.8 and 1 (Moe and
Moe, 2005). The spacecraft surface temperature will be assumed
R
Cl dA constant at Tw = 300 K.
CL = (26) Note that the drag and lift coefficients are dependent on the
Aref
atmospheric parameters through the vre /V∞ and s parameters.
As the atmosphere has temporal and spatial variability (vertical
γ = cos (ϕ) (27)
but also horizontal) the force coefficients will in general be
variable during an orbit. These changes in the force coefficients
l = sin (ϕ) (28) are small and can be safely ignored given that the variability
of the atmosphere (changes in atmospheric density and relative
1
G= (29) flow direction and magnitude) is orders of magnitude larger.
2s2 Additionally, although the Sentman model can provide the lift
coefficient CL , it is, in general, an order of magnitude smaller than
1 −γ 2 s2
P= e (30) the drag coefficient CD and thus it will be neglected in this study
s (Doornbos, 2011).
Q=1+G (31) • Asm.10: The changes of vre /V∞ and s during an orbit
are negligible when compared to the atmospheric
Z = 1 + erf (γ s) (32) density variability.

The ϕ term denotes the angle between the flow and the surface
4.4. Aerodynamic Properties of a Sphere
local normal vector (0◦ when the surface is normal to the flow
Equation (34) can be used to compute the drag coefficient CD of
and 90◦ when it is parallel), vre the most probable velocity of
a sphere. The reference area is set as the cross section area of the
the reemitted gas particles, V∞ the relative bulk velocity between
sphere Aref = πR2 . In Equation (34), θsc and φsc are the azimuthal
the spacecraft and the incident gas particles (the same one as in
and polar spherical coordinate angles.
Equations 20–22), Aref an arbitrary reference area (usually the
cross section area of the spacecraft), s the ratio between V∞ and R R π R 2π
the most probable thermal velocity of the gas vth (s = Vv ∞ ), and Cd (ϕ) dA 0 0 Cd (θsc , ψsc ) sin φsc dθsc dφsc
th CD = = (34)
erf (x) denotes the error function. πR2 π
According to Koppenwallner (2009) the vre /V∞ ratio can be
written as in Equation (33), with Rg denoting the gas constant Equation (35) computes the angle between the flow and the local
and Tw the temperature of the surface (wall). normal vector ϕ (required by the Sentman model) using the polar
s  spherical coordinate angles.
4Rg Tw
 
vre 1
= 1 + σa 2
− 1 (33)
V∞ 2 V∞ cos ϕ = cos (π/2 − φsc ) cos θsc (35)

Frontiers in Robotics and AI | www.frontiersin.org 7 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

Due to the sphere’s symmetry, the drag coefficient is constant The location of the spacecraft CoP in body axes pCoP can then be
regardless of the orientation of the sphere with respect to the written as in Equation (38).
flow (greatly simplifying the analysis). By using Equation (34) a    
sphere’s drag coefficient is found to be around CD ≈ 2.1. Figure 4 −dCoM0 dCoP
clearly shows how the drag coefficient of a sphere changes with pCoP = 0  + RBF  0  (38)
altitude, solar activity and energy accommodation coefficient, 0 0
making this magnitude variable and uncertain.
The orientation of the body reference frame B0 with respect Due to the sphere’s symmetry, the relative flow îflow , the CoP,
to the orbital frame will be denoted by the common roll φ, pitch and the sphere’s geometric center are aligned. Therefore, there
θ and yaw ψ Euler angles. When roll, pitch and yaw are 0 the is no torque with respect to the sphere’s geometric center. The
body frame is aligned with the orbit frame. The relative flow aerodynamic torque with respect the host vehicle CoM is then
direction will be defined with its own reference frame where only a function of dCoM0 as shown in Equation (39). It can then
the flow direction will be in −îflow axis. The orientation of this be assumed that torque-wise, the effective CoP is located at the
flow reference frame will be denoted by a flow pitch θflow and sphere’s geometric center.
flow yaw ψflow .
Let ROB denote the rotation matrix from body to orbital and
 
−dCoM0
ROF the rotation matrix from flow to orbital and thus RBF = τ aero = pCoP × F aero = 0  × F aero (39)
RTOB ROF is the rotation from the flow to the body reference frame. 0
With these definitions the aerodynamic force in body axes is
defined by Equation (36), with D denoting the aerodynamic drag. Although in this work a spherically shaped spacecraft has been
used, an analogous analysis can be conducted for spacecraft with

−D
 more complex shapes.
F aero = RBF  0  (36) It may be useful when devising the controllers to simplify these
0 aerodynamic force and torque equations. If the control method
fulfills its goal the spacecraft attitude will be in close vicinity
of its target attitude φ ≈ 0, θ ≈ 0, ψ ≈ 0 (small angles
The CoP of the sphere is aligned with the direction of the flow approximation). Additionally, the atmospheric co-rotation and
îflow (see Figure 5). The location of the center of pressure from wind do not cause the relative flow to have large deviations with
the sphere’s center can be computed using (Equation 37). Again respect to the inertial velocities (see Figure 3B) making θflow
this magnitude is slightly dependent on the altitude but it will be and ψflow also small. Under these assumptions, the Euler angles
assumed as constant (it will eventually be shown that the location of the spacecraft with respect to the flow (rotation represented
of the real CoP of the sphere is not relevant). by RBF ) can be approximated using φ ′ = −φ, θ ′ = θflow −
θ and ψ ′ = ψflow − ψ (which will also be small angles)
and the aerodynamic forces in body axes can be subsequently
approximated by Equation (40).
R
dCoP Cd (ϕ) xdA
=
R πR2 CD    
R π R 2π cos (θflow − θ ) cos (ψflow − ψ) −1
Cd (θsc , φsc ) sin2 φsc cos θsc dθsc dφsc F aero ≈ −D  cos (θflow − θ ) sin (ψflow − ψ)  ≈ D  −ψ ′ 
= 0 0 ≈ 0.66
πCD − sin (θflow − θ ) cos (ψflow − ψ) θ′
(37) (40)

FIGURE 4 | Variation of a sphere drag coefficient with altitude, solar activity (A) and energy accommodation coefficient (B).

Frontiers in Robotics and AI | www.frontiersin.org 8 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

A simplified expression for the aerodynamic torque can also be aerodynamic torques will provide a restoring torque if the sphere
obtained using the same small angle approximation as shown in center is behind the host vehicle CoM dMp > 0, making the
Equation (41). spacecraft oscillate around this equilibrium point (marginally
  stable). If the host vehicle CoM is leading the center of the sphere
0 dMp < 0 the system is unstable.
τ aero ≈ D  θ ′ dCoM0  (41) In a marginally stable configuration dMP > 0 the natural
ψ ′ dCoM0 frequency of the oscillation can be approximated (using the small
angles approximation) with Equation (42).
From Equation (41) it can be clearly seen that the equilibrium r
attitude is that attitude where the flow, the host vehicle CoM, DdCoM0
and the sphere’s geometric center are aligned (θ ′ = ψ ′ = 0). ωn = (42)
I
When there is a misalignment of this equilibrium attitude, the
The natural frequency will, in general, be small and thus it can
be normalized with the orbital mean motion to make it easier
to read. Figure 6 shows the natural frequency (normalized with
the orbit period) for two different sphere sizes (R = 10 cm and
R = 25 cm) with different CoM to CoP distances, at different
altitudes and using the numerical parameters shown in Table 1.
As the magnitude of the aerodynamic force is proportional
to the area ∼ R2 , larger spacecraft with equivalent natural
frequencies will exhibit smaller perturbations as will have larger
inertias ∼ R5 . Another expected result is that as the altitude
increases and the aerodynamic disturbance weakens, the natural
frequency also decreases. Thus, it is readily apparent that the
aerodynamic disturbances will be more important for small
spacecraft at low altitudes. The simulations that have been
conducted have then been focused on small spacecraft examples.
It is important to note that by definition, the host vehicle
CoM is displaced with respect to the sphere center only along
the B0 î direction. This condition has been imposed to simplify
the analysis but it is expected to be met by VLEO spacecraft. In a
generic case, the CoM can be displaced in any direction and then
a secular aerodynamic torque will appear when the spacecraft
is at the target attitude (ignoring the direction variability of the
relative flow direction). Therefore, it is highly desirable to have
FIGURE 5 | Location of the sphere’s Center-of-Pressure.
the CoM and effective CoP (center of the sphere) aligned with
the B0 î axis in order to avoid these secular torques. As it is

FIGURE 6 | Natural frequency (normalized with the orbit period) for different CoM and altitudes for R = 10 cm (A) and R = 25 cm (B).

Frontiers in Robotics and AI | www.frontiersin.org 9 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

TABLE 1 | Numerical parameters. • Reduced Model Controller Asm.2: Shifting mass is a point mass.
• Reduced Model Controller Asm.3: The mass and inertia
Parameter Value
properties of the host vehicle and of the shifting mass
κ 0.1 are known.
ρs 500 kg m2 • Reduced Model Controller Asm.4: The relative position,
CD 2.2 velocity and acceleration of the shifting mass are known.
σC D 0.22 (10 % of the nominal CD = 2.2)
Solar activity indices Moderate activity as in ISO 14222 (2013).
Under these assumptions, the equation of motion in Equation
(19) can be further simplified, yielding Equation (43). The
position, velocity and acceleration of the unique shifting mass
expected that VLEO spacecraft designers will take this issue into with respect to the body axes is denoted by x, y and the shifting
consideration it can be safely assumed that the CoM and CoP, mass velocity with respect to the body reference frame B0
for the target attitude, will be reasonably aligned. Any residual by ẋ′ , ẏ′ .
misalignment can be corrected by a bias in the position of the
Jz + µ x2 + y2 ψ̈ + µ 2 xẋ′ + yẏ′ ψ̇ + xÿ′ − yẍ′

 

shifting masses (resulting in a reduction of their shifting range
µ 
and control authority) and thus the initial assumption of host

= τz + Fx y − Fy x (43)
vehicle CoM displaced only along the B0 î direction can be M0
recovered (without loss of generality).
The aerodynamic disturbances have low frequencies (similar
The mission of the shifting masses is then to stabilize
to the orbit frequency) and so it is expected that the motion
the spacecraft in the presence of these aerodynamic attitude
of the shifting mass will be also slow (small velocities and
disturbance eliminating the need to use other actuators for this
accelerations), thus limiting the dynamic effects of the shifting
purpose (thus potentially delaying saturation and saving power
mass. Additionally, as the shifting mass m is small compared to
and mass).
the host vehicle mass µ ≪ 0, the dynamic effects of the shifting
The assumption that the host vehicle CoM is displaced
mass will be further reduced and they can therefore be safely
with respect to the sphere center roughly along the B0 î
neglected during the controller design.
direction represents one of the limitations of the proposed
method. The shifting masses are only able to reject the • Reduced Model Controller Asm.5: Shifting mass velocities and
aerodynamic disturbances for the limited set of attitudes where accelerations have negligible effects on the dynamics.
this assumptions holds. For arbitrary attitudes, the aerodynamic
As the shifting range is also small the change on the system’s
torques may be too strong to be compensated by the shifting
inertia is also small and thus the system’s inertia will be
masses. However, in the context of VLEO, these strong secular
considered as constant (using the initial shifting masses position
torques may be a burden for other attitude control methods as
x0 and y0 ) during the controller design. These assumptions
well (e.g., rapidly saturating reaction wheels) and thus holding
further simplify the equations of motion to Equation (44). It also
attitudes far from the aerodynamic equilibrium points is an
has to be noted that only aerodynamic forces and torques will
intrinsic challenge for spacecraft operating in VLEO.
be considered.
µ 
5. REJECTION CAPABILITIES UNDER A Jz + µ x02 + y02 ψ̈ = τz +

 
Fx y − Fy x (44)
M0
REDUCED MODEL WITH ONE
ROTATIONAL DEGREE-OF-FREEDOM
• Reduced Model Controller Asm.6: Constant system inertia.
A potential application of the shifting masses is disturbance
rejection. For simplicity, it is worth to start the rejection The shifting masses modulate τ̂ aero by actively changing the
capability analysis with a reduced model that only considers location of the system CoM. Using the aerodynamic properties
one rotational degree-of-freedom. This analysis provides insight of a sphere and using the small angles approximation, Equation
into the rejection capabilities and the shifting mass system (45) can be obtained.
requirements with respect to the system’s parameters. It is also  
µ 
Jz + µ x02 + y02 ψ̈ = D ψ ′ dCoM0 + −y + ψ ′ x

 
of particular interest to explore how the spacecraft size and
operating altitude drive the required shifting mass and range in M0
(45)
order to meet pre-specified performance requirements.
The yaw ψ rotation has been selected for this one-dimensional
analysis as the co-rotation and predominant wind act on this • Reduced Model Controller Asm.7: The system remains at all
particular axis. Additionally, a single shifting point mass will be times close to its target attitude (small angles approximation).
used and the controller will be based on linearized dynamics. The
It is immediately clear from Equation (45) that to generate a
goal of the shifting mass is then to stabilize the spacecraft around
control torque it is much more effective for the mass to move
ψ = 0 and reject the disturbance induced by ψflow .
perpendicular to the relative flow (in this case y) than parallel to
• Reduced Model Controller Asm.1: Single rotational it (along x). So in order to limit the system complexity, it will be
degree-of-freedom (yaw axis). assumed that the shifting mass only moves in y (perpendicular

Frontiers in Robotics and AI | www.frontiersin.org 10 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

to the flow if ψ ′ is small). Shifting the mass only along one To explore the design space and the system response it will
direction reduces the volume and the complexity of the shifting be assumed that the PID controller is tuned so that the closed
mass system while maximizing its effectiveness. It is understood loop system has a specific bandwidth and phase margin. In
that (Equation 45) has been simplified for small angles and thus these fixed controller conditions a Montecarlo simulation can
shifting the mass along y will only be perpendicular to the flow be performed to extract the required shifting mass range for a
direction only for ψ = ψflow = 0, if there is a large misalignment given spacecraft size and the uncertain aerodynamic properties
the y shifting mass will start to loose efficacy. and environmental conditions.
Although the controller is build upon a linearized model (see
• Reduced Model Controller Asm.8: The shifting mass only
all Reduced Model Controller Asm.), the numerical simulations
moves along the ĵ axis.
use the full dynamic equations and the high-fidelity environment
Another important consideration that it is apparent from models (only using the generic Asm.). To emulate the uncertainty
Equation (45) is that the maximum torque provided by the on the aerodynamic properties, the actual drag coefficient used in
shifting mass is τmax = ±D M0m+m ymax . It is clear that the mass the numerical simulation differs from the one used to design the
of the shifting mass and the available shifting range are the two controller. Although the CoP is considered known, given that a
variables at the designer disposal to regulate the control authority spherical shape is used, the uncertainty on the drag coefficient
of the system. can also be used to emulate an uncertainty in the CoP position.
The atmospheric density and the magnitude and direction of Sun-synchronous circular orbits with 10:30 am mean LTAN have
the flow are inherently unknown to the controller. An estimated been used.
density, purely based on the altitude (no horizontal variability) The Montecarlo simulations are initialized with ideal stable
will be used by the controller. Additionally, the controller will attitudes ψ = 0 and ψ̇ = 0 and thus emulate steady state
assume that the incident flow matches the inertial velocity conditions. Each Montecarlo run simulates 4 consecutive orbits
magnitude and direction. and 100 simulations are used to extract the statistics (with error
bars denoting the 95% confidence interval).
• Reduced Model Controller Asm.9: Constant
Figure 7 shows the maximum shifting range and attitude error
atmospheric density.
(3σ values) for a 10 cm radius spherical satellite for different
• Reduced Model Controller Asm.10: Relative flow velocity
mass fractions of the shifting mass m/M0 and for a CoM leading
matches the spacecraft’s inertial velocity.
the CoP by 3% of the sphere radius R. Figure 8 shows how the
Under these conditions, the system equations can be written required shifting range and attitude error change for different
as Equation (46) which corresponds to the transfer function CoP to CoM distances dCoM0 and with a fixed shifting mass
written in Equation (47). This represents a simple second fraction kept at 3% of the host vehicle mass M0 .
order Single Input Single Output system and a Proportional- The system bandwidth in the controllers used to generate
Integral-Derivative (PID) controller can be easily designed and (Figures 7,8) has been kept at four times the natural frequency
implemented to reject the aerodynamic disturbances while of the system 4ωn and the phase margin set to 30 deg. This allows
keeping the spacecraft stable at ψ = 0. a comparison of the system performance even if the altitude or
    spacecraft size are changed.
 2 2

 µ µ For the 10 cm radius case it is quite clear that the proposed
Jz + µ x0 + y0 ψ̈ = D −ψ x0 + dCoM0 − y
M0 M0 method is able to reject the aerodynamic disturbances and
(46) maintain a reasonably stable attitude (with respect to 10 cm sized
spacecraft standards Polat et al., 2016) with mass fraction and
ψ (s) b
T (s) = = ′2 (47) shifting range requirements compatible with the spacecraft mass
y (s) J s +k and dimension constrains (considering that realistic uncertainty
in the parameters has been taken into account). Note how the
J ′ = Jz + µ x02 + y02 proposed method is able to stabilize the spacecraft even if the


(48)
CoM is behind the CoP (unstable configuration). As expected the
  required mass fraction and required shifting range decrease as
µ
k=D x0 + dCoM0 (49) the CoP gets closer to the CoM. It is also worth pointing that the
M0 unstable configuration dp < 0 requires higher shifting range than
their stable counterparts.
µ The attitude error, which is constant in Figure 7 due to the
b = −D (50)
M0 constant bandwidth employed, can be decreased if the bandwidth
of the close looped system is increased. The required shifting
This controller also carries the underlying assumption that
range can be decreased by decreasing the phase margin. But
the shifting mass can instantaneously move, without lag,
both measures have limits. By decreasing the phase margin
from one position to another one. This will be relaxed in
the controller is less robust and increasing the bandwidth
subsequent sections.
increases the gains which imposes more strict requirements
• Reduced Model Controller Asm.11: Shifting mass movement on the sensors and actuators. The PID gains for the 10 cm
has infinite bandwidth. radius case are shown in Figure 9 when the angular and

Frontiers in Robotics and AI | www.frontiersin.org 11 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

FIGURE 7 | Required 3σ shifting range (A) and 3σ attitude error (B) with respect to the shifting mass fraction for a 10 cm radius spherical spacecraft.

FIGURE 8 | Required 3σ shifting range (A) and 3σ attitude error (B) with respect to the CoP to CoM distance dCoM0 for a 10 cm radius spherical spacecraft.

angular velocity errors are provided in rad and rad/s. As Lower natural frequencies do reduce the proportional and
expected, the gains increase as the system becomes more difficult derivative gains and thus spacecraft of bigger size or operating at
to control, that is with increasing dCoM0 or reducing mass higher altitudes have a wider margin to increase their controller
fraction of the shifting mass. It is important to note that bandwidth and reduce the attitude error whilst maintaining
the high derivative gain may impose certain requirements on reasonable gains.
the attitude velocity estimates that the attitude determination The tuning employed in this examples appears to
subsystem needs to provide to the controller (specially in terms of give satisfactory results with the selected parameters and
noise levels). uncertainties. These examples illustrate the general trends
It is also interesting to note that the relative shifting range and provide performance estimates that can be later used as
increases slightly with altitude or spacecraft size as a relative initial guesses.
larger shifting is required to generate the same acceleration at
higher altitude or for spacecraft with larger inertias. As in this
analysis the controller bandwidth is kept constant relative to the 6. THREE ROTATIONAL
natural frequency, larger spacecraft display larger attitude errors DEGREES-OF-FREEDOM CASE
(as their natural frequency is significantly lower). For a 25 cm
radius spherical satellite the attitude error and required shifting The previous analysis has been conducted using a reduced model
range for different altitudes, mass fractions, and CoP to CoM and only considering a single rotational degree-of-freedom. That
distances is shown in Figures 10, 11. analysis is useful as provides generic results and shows the trends

Frontiers in Robotics and AI | www.frontiersin.org 12 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

FIGURE 9 | PID gains (A) for different mass fractions and (B) for different CoP to CoM distances dCoM0 for a 10 cm radius spherical spacecraft.

FIGURE 10 | Required shifting range (A) and attitude error (B) with respect to the shifting mass fraction for a 25 cm radius spherical spacecraft.

FIGURE 11 | Required shifting range (A) and attitude error (B) with respect to the CoP to CoM distance for a 25 cm radius spherical spacecraft.

Frontiers in Robotics and AI | www.frontiersin.org 13 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

when the different parameters are varied. In this section two and their mass fractions are small. When building this controller
simple controllers for the three rotational degree-of-freedom case many of the assumptions made for the PID one rotational degree-
will be presented to demonstrate the applicability of the proposed of-freedom controller are carried over.
method to a multidimensional case.
• LQR Controller Asm.1: Linearized dynamics.
6.1. Shifting Masses Driver • LQR Controller Asm.2: Shifting masses are point masses.
For the PID controller it was assumed that the shifting masses • LQR Controller Asm.3: The mass and inertia properties of the
motion had infinite bandwidth. As the movement of the masses host vehicle and of the shifting masses are known.
was slow that did not pose any problem. In reality the controller • LQR Controller Asm.4: The relative position, velocity and
will specify the position of shifting masses and then the accelerations of the shifting masses are known.
masses will need to move to those positions using a limited • LQR Controller Asm.5: Shifting masses velocities and
acceleration and velocity. For the three rotational degrees-of- accelerations have negligible effects on the dynamics.
freedom controllers an underlying PID controller will control the • LQR Controller Asm.6: Constant system inertia.
motion of the shifting masses to their desired location. • LQR Controller Asm.7: The system remains at all times close to
Although the underlying shifting masses actuators may be its target attitude (small angles approximation).
capable of rapid motion and aggressive acceleration, these • LQR Controller Asm.8: Constant atmospheric density.
magnitudes may need to be bounded in order to limit the • LQR Controller Asm.9: Relative flow velocity matches the
dynamic effects of the shifting masses motion (which are spacecraft’s inertial velocity.
currently ignored by the controllers and thus could potentially When the spacecraft is in the vicinity of its equilibrium point
degrade the controller’s performance). This is specially relevant the masses that shift perpendicular to the relative flow provide
when the shifting masses controller may request abrupt shifting the maximum efficacy. As the goal is to keep the spacecraft
masses position changes (e.g., when alternating saturation stable then using only two shifting masses (with mass m1 and
positions are requested). m2 ), respectively moving along the B0 pitch ĵ and yaw k̂ axes,
Equation (43) can be used to select these velocity and maximizes the available torque while minimizing the system’s
acceleration limits. For example these limits can be selected so complexity and the required volume.
that the effects due to the shifting masses velocity or acceleration
is no greater than a certain fraction η of the aerodynamic torque • LQR Controller Asm.10: Two masses moving along the ĵ and
created by that shifting mass, as shown in Equations (51, 52). k̂ axes.

D As the shifting masses are unable to provide any control torque

′
ẏmax <η (51) parallel to the flow then an additional third actuator is required.
2ψ̇max M0 In this analysis an ideal actuator acting on roll î providing τroll
has been assumed.
′ Dymax • LQR Controller Asm.11. Ideal roll actuator to augment the
ÿmax <η (52)
M0 x0 otherwise underactuated system.
These proposed limits are mainly dominated by the host vehicle In the typical case where the system is in the vicinity of the
mass M0 and by the atmospheric density (through the drag equilibrium point, this configuration and the LQR type controller
D). As the host vehicle mass is proportional to the radius of is well suited. For a detumbling case, where the system does
the spacecraft these limits can also be set with respect of the not remain in the linear region, an LQR controller may be
vehicle size and in this analysis a simple shifting mass acceleration less efficient than other non-linear controllers. It is also worth
and velocity limit of R/10 m/s2 and R m/s respectively has mentioning that, as the aerodynamic force is the dominant
been employed. perturbation and only acts perpendicular to the flow direction
Another option which would regulate itself would be to (nominally along pitch and yaw), it is expected that the roll
include the dynamic effects of the shifting masses motions actuator will not be required to provide significant torques when
into the controller. The controller would then avoid sudden the system is in the stable around the target equilibrium attitude.
accelerations or try to compensate the dynamic effects caused by As in all LQR based controllers, the weight matrices have to be
the shifting mass motion. carefully adjusted to obtain the desired combination of attitude
error, shifting mass range and sufficiently small gains (so that the
6.2. Linear Quadratic Regulator Approach requirements on the sensors can be met). Figures 12, 13 show a
Given the good performance of the PID controller on a steady state example of the evolution of the attitude, the shifting
one rotational degree-of-freedom case, it seems that a Linear masses displacement and roll torque when an LQR controller is
Quadratic Regulator (LQR) based controller may also have a used and when the gravity gradient perturbation is also included.
good performance for an attitude hold scenario (steady state) in A 25 cm spherical satellite operating at a 300 km Sun-
a three rotational degrees-of-freedom case. When the system is Synchronous and with 6% of the mass allocated to the shifting
linearized around the equilibrium point all the non-linearities masses (3% for each mass) and an estimated distance between the
from the presence of the shifting masses disappear or are CoP and the host vehicle CoM of 3% of the radius have been used.
neglected as it is assumed that the shifting masses move slowly In this case the spacecraft is operating in the unstable orientation

Frontiers in Robotics and AI | www.frontiersin.org 14 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

the angular velocity error (for a stable attitude with respect to the
orbit frame a pitch that matches the orbital motion needs to be
included), J is the inertia matrix and the K p and K d diagonal
matrices being the proportional and derivative gains in each
axis respectively.
Some guidance to select the gains can be obtained by
assuming small angles, a single degree-of-freedom, and that the
shifting mass motion dynamics effects are negligible. In that
case, the system reduces to a second order system and thus the
proportional K p and derivative gains K d for each axis (kp and kd )
can be related to the desired closed loop natural frequency ωn and
a damping ratio ξ as shown in Equation (54) (Wie et al., 1989).

kp = 2ωn2 kd = 2ξ ωn (54)

The quaternion feedback is particularly suited to be used when

large attitude misalignments are present and thus it will be
FIGURE 12 | Roll, pitch and yaw using the LQR controller. employed here for a detumbling scenario. If it is assumed that the
system is composed by three shifting masses each moving along
(CoP leading the CoM). Uncertainty in the drag coefficient CD the body B0 roll î, pitch ĵ and yaw k̂, the aerodynamic torque
has also been included. Although the controller is based on the provided by the shifting masses can be written as follows (see
linearized equations of motion, in the simulator the full equations Equation 17).
of motion have been used. As seen in Figures 12, 13 the controller  
is able to maintain the system stable with small shifting masses m 1 r1
F aero
and shifting range. τ sm = ×  m 2 r2  (55)
M + m1 + m2 + m3 m r 3 3
6.3. Quaternion Feedback With Partial
Feedback Linearization With m1 , m2 , and m3 denoting the shifting masses and r1 , r2 and
A more general approach that does not rely on linearization r3 their shifting ranges along î, ĵ and k̂ respectively. Note how
uses the well known quaternion feedback with a partial feedback Equation (55) is the last term of Equations (16, 17).
linearization approach (Wie and Barba, 1985; Wie et al., 1989). As the controller has no information about the actual direction
The estimated aerodynamic torque τaero is also used to help and magnitude of the aerodynamic force, an estimate F aero needs
with the feedback linearization but the terms related the shifting to be used in the steering logic. A relative flow matching the
masses motion are left out as they are assumed to be negligible. inertial velocity and a mean atmospheric density are used to
Again, many of the assumptions are carried over. obtain this estimate. It is clear from Equation (55) that the
shifting masses aerodynamic torque is perpendicular to the F aero
• Quat. Feedback Controller Asm.1: Shifting masses are and thus other actuators should provide the required torque that
point masses. is parallel to F aero .
• Quat. Feedback Controller Asm.2: The mass and inertia The shifting masses position to achieve the requested
properties of the host vehicle and of the shifting masses perpendicular torque can be obtained using (Equation 56). Note
are known. that (Equation 56) is equivalent to the expression used to drive
• Quat. Feedback Controller Asm.3: The relative position, the magnetic torquers (replacing the magnetic field with the
velocity and accelerations of the shifting masses are known. aerodynamic force F aero and the magnetic moment with the
• Quat. Feedback Controller Asm.4: Shifting masses velocities mr product).
and accelerations have negligible effects on the dynamics.
• Quat. Feedback Controller Asm.5: Constant system inertia.
 
m 1 r1 τ × F aero (M + m1 + m2 + m3 )
• Quat. Feedback Controller Asm.6: Constant  m2 r2  = req (56)
atmospheric density. m 3 r3 F aero · F aero
• Quat. Feedback Controller Asm.7: Relative flow velocity
matches the spacecraft’s inertial velocity. The resulting shifting masses provided torque τ sm is just
Under this control law the requested torque can be written as in the component perpendicular to the flow direction and then
Equation (53). additional actuators need to provide the parallel torque. To keep
it consistent with past numerical examples only two shifting
τ req = −K p Jqe − K d Jωe + ω0 × Jω0 − τ aero (53) masses shifting along the body’s B0 pitch ĵ or yaw k̂ axes will
be used (thus r1 = 0). These couple of shifting masses will
In Equation (53) qe denotes the vector elements of the error be complemented by a single ideal actuator acting along the B0
quaternion (Wie and Barba, 1985; Wie et al., 1989), ωe denotes roll î axis.

Frontiers in Robotics and AI | www.frontiersin.org 15 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

FIGURE 13 | Shifting (A) and roll torque (B) required for an LQR controller.

7. PRACTICAL IMPLEMENTATION ON THE

SHIFT-MASS SAT 3U CUBESAT
Small spacecraft are more sensitive to aerodynamic disturbances
due to their high area to inertia ratio. A CubeSat operating
at low altitude is thus a good first candidate to implement
the proposed aerodynamic disturbance rejection method. This
implementation exercise using a CubeSat, which are highly
constrained platforms in terms of mass and volume, also serves
as a practical feasibility check of the whole concept.
A preliminary design of the “Shift-Mass Sat” 3U CubeSat with
three orthogonal shifting masses is shown in Figure 15. All the
components, subsystems and shifting masses, are COTS to ensure
their commercial availability. More detailed information on this
design can be found in Polat (2016).
The three 150 g shifting masses approximately take 75%
of a 1U volume and have a 70 mm useful travel range.
Magnetic torquers augment the shifting masses and complete the
FIGURE 14 | 3σ maximum stabilization time for a 25 cm radius spacecraft in a
400 km altitude orbit. actuator set.
The performance of this design, paired with an LQR
controller, has been evaluated using a numerical simulation.
The CoM position, mass, inertia, shifting mass and travel
• Quat. Feedback Controller Asm.8: Two shifting masses moving
range, magnetic dipole moment of the magnetic torquers, and
along the ĵ and k̂ axes .
aerodynamic properties used for the simulation have been
• Quat. Feedback Controller Asm.9. Ideal roll actuator to
derived from the prototype design. It is worth pointing that for
partially provide the required torque parallel to F aero .
this particular design, the shifting masses have a 1.82 mm control
This configuration has been used to stabilize a 25 cm spacecraft authority on the combined system’s CoM position to modulate
with a random initial attitude and an initial angular velocity the aerodynamic torque direction and magnitude.
in a random orientation with a magnitude between ±0.5 deg/s The LQR controller is used for detumbling and to keep the
(chosen with a random uniform distribution). Figure 14 shows spacecraft stable. A gain scheduling scheme, where less aggressive
the mean and the 3σ maximum stabilization time obtained by gains are employed during the detumbling phase is employed.
a 25 sample MonteCarlo simulation. The control law gains have The initial angular velocity of the CubeSat is chosen as 0.01
been set according to Equation (54) with a bandwidth of twice rad/s in all axis and the orbit altitude is set at 300 km. The
the spacecraft natural frequency and ξ = 0.7. The shifting masses shifting masses movement and Euler angles of one of these
represent 6% of the host vehicle mass (3% each shifting mass) and simulations are presented in Figure 16 showing the feasibility of
the CoM leads the CoP by 3% of the spacecraft radius. the proposed method.

Frontiers in Robotics and AI | www.frontiersin.org 16 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

FIGURE 15 | Prototype 3U CubeSat design with shifting masses (A) and detail of the three orthogonal shifting masses (B).

FIGURE 16 | Euler angles (A) and shifting masses positions (B).

It is also worth noting that in Figure 16B the shifting masses Commercial-Off-the-Shelf components and an Linear Quadratic
exhibit a bias when the spacecraft is stabilized. This is because the Regulator controller demonstrates its technological feasibility.
CubeSat center of mass is not completely centered around the Further research could be directed to develop other types
roll axis. of controllers, specially non-linear controllers, to drive the
shifting masses in order to increasing the performance of the
8. CONCLUSIONS system. However, rigorously proving their stability can be a
challenging endeavor.
In conclusion, using a set of shifting masses that shift the
spacecraft’s center-of-mass is a viable method to reject the AUTHOR CONTRIBUTIONS
aerodynamic disturbances present at Very Low Earth Orbit.
Despite the highly non-linear dynamics of a spacecraft with JV-L and HP did most of the research. JV-L wrote the manuscript.
internal moving parts simple controllers based on the linearized MR supervised the research.
equations of motion suffice to keep the spacecraft stable. The
requirements imposed on the attitude determination subsystem ACKNOWLEDGMENTS
and the shifting masses (shifting range and mass fraction) are
well within practical limits. Achieving stabilization from arbitrary Preliminary results of this research effort were presented at the
initial attitude and small angular velocities is also possible. 26th AAS/AIAA Space Flight Mechanics Meeting in Napa, CA,
A prototype implementation on a 3U CubeSat only using February 14-18, 2016 (Virgili-Llop et al., 2016).

Frontiers in Robotics and AI | www.frontiersin.org 17 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

REFERENCES Koppenwallner, G. (2009). Energy Accommodation Coefficient and Momentum

Transfer Modeling. HTG–TN–08–11, HTG (Katlenburg Lindau).
Ahn, Y. T. (2012). Attitude Dynamics and Control of a Spacecraft Using Shifting Kumar, K. (2010). “Attitude control of miniature satellites using movable masses,”
Mass Distribution. Ph.D. thesis, The Pennsylvania State University. in Ryerson University SpaceOps 2010 Conference (Huntsville, AL: American
Atkins, B. M., and Henderson, T. A. (2012). Under-actuated moving mass attitude Institute of Aeronautics and Astronautics).
control for a 3u cubesat mission. Adv. Astronautical Sci. 143, 2083–2094. Kumar, R. R., Mazanek, D. D., and Heck, M. L. (1995). Simulation and shuttle
Bevilacqua, R., and Romano, M. (2008). Rendezvous maneuvers of multiple hitchhiker validation of passive satellite aerostabilization. J. Spacecr. Rockets 32,
spacecraft using differential drag under J2 perturbation. J. Guid. Control Dyn. 806–811. doi: 10.2514/3.26688
31, 1595–1607. doi: 10.2514/1.36362 Kumar, R. R., Mazanek, D. D., and Heck, M. L. (1996). Parametric and classical
Bird, G. A. (1994). Molecular Gas Dynamics and the Direct Simulation of Gas Flows. resonance in passive satellite aerostabilization. J. Spacecr. Rockets 33, 228–234.
Oxford, UK: Oxford Science Publications. doi: 10.2514/3.26745
Bowman, D. S., and Lewis, M. J. (2002). Minimum drag power-law shapes for Kunciw, B., and Kaplan, M. (1976). Optimal space station detumbling by
rarefied flow. AIAA J. 40, 1013–1015. doi: 10.2514/2.1743 internal mass motion. Automatica 12, 417–425. doi: 10.1016/0005-1098(76)
Challinor, R. (1968). The apparent rotation of the upper atmosphere. Planet. Space 90003-0
Sci. 16, 557–566. doi: 10.1016/0032-0633(68)90097-4 Larsen, M. F., and Fesen, C. G. (2009). Accuracy issues of the existing
Chesi, S., Gong, Q., and Romano, M. (2017). Aerodynamic three-axis attitude thermospheric wind models: can we rely on them in seeking
stabilization of a spacecraft by center-of-mass shifting. J. Guid. Control Dyn. solutions to wind-driven problems? Ann. Geophys. 27, 2277–2284.
40, 1613–1626. doi: 10.2514/1.G002460 doi: 10.5194/angeo-27-2277-2009
Doornbos, E. (2011). Thermospheric Density and Wind Determination From Moe, K., and Moe, M. M. (2005). Gas-surface interactions and satellite drag
Satellite Dynamics. Ph.D. thesis, Technische Universiteit Delft. coefficients. Planet. Space Sci. 53, 793–801. doi: 10.1016/j.pss.2005.03.005
Drinkwater, M. R., Haagmans, R., Muzi, D., Popescu, A., Floberghagen, R., Kern, Moe, K., and Moe, M. M. (2010). “Gas-surface interactions in Low-Earth Orbit,” in
M., et al. (2007). “The GOCE gravity mission: ESA’s first core earth explorer,” in 27th International Symposium on Rarefied Gas Dynamics, 2010 (Pacific Grove,
Proceedings of the 3rd International GOCE User Workshop, 6-8 November, 2006 CA: American Institute of Physics).
(Frascati: ESA), 1–8. Moe, K., Moe, M. M., and D, S. (1998). Improved satellite drag coefficient
Drob, D. P., Emmert, J. T., Crowley, G., Picone, J. M., Shepherd, G. G., Skinner, W., calculations from orbital measurements of energy accommodation. J. Spacecr.
et al. (2008). An empirical model of the Earth’s horizontal wind fields: HWM07. Rockets 35, 266–272. doi: 10.2514/2.3350
J. Geophys. Res. Space Phys. 113. doi: 10.1029/2008JA013668 Pardini, C., Moe, K., and Anselmo, L. (2012). Thermospheric density model
ECSS Secretariat (2008). ECSS Space Engineering-Space Environment. Technical biases at the 23rd sunspot maximum. Planet. Space Sci. 67, 130–146.
Report ECSS-E-ST-10-04C, ESA. doi: 10.1016/j.pss.2012.03.004
Edwards, T. L., and Kaplan, M. H. (1974). Automatic spacecraft detumbling by Picone, J. M., Hedin, A. E., Drob, D. P., and Aikin, A. C. (2002).
internal mass motion. AIAA J. 12, 496–502. doi: 10.2514/3.49275 NRLMSISE-00 empirical model of the atmosphere: statistical comparisons
Fortescue, P. W., and Stark, J. P. W. (1995). Spacecraft Systems Engineering, 4th and scientific issues. J. Geophys. Res. Space Phys. 107, SIA 15-1–SIA 15-16.
Edn. Chichester, UK: Wiley and Sons, Inc. doi: 10.1029/2002JA009430
Gargasz, M. L. (2007). Optimal Spacecraft Attitude Control Using Aerodynamic Polat, H. C. (2016). Prototype Design and Mission Analysis for a Small Satellite
Torques. Master’s thesis, Air Force Institute of Technology. Exploiting Environmental Disturbances for Attitude Stabilization. Master’s
Gregory, J. C., and Peters, P. N. (1987). “A measurement of the angular distribution thesis, Naval Postgraduate School.
of 5 ev atomic oxygen scattered off a solid surface in Earth orbit,” in Proceedings Polat, H. C., Virgili-Llop, J., and Romano, M. (2016). Survey, statistical analysis
of the 15th International Symposium on Rarefied Gas Dynamics, Vol. 1, eds V. and classification of launched cubesat missions with emphasis on the attitude
Boffi and C. Cercignani (Grado), 644–656. control method. J. Small Satellites 5, 513–530.
Grubin, C. (1962). Dynamics of a vehicle containing moving parts. J. Appl. Mech. Prieto, D. M., Graziano, B. P., and Roberts, P. C. (2014). Spacecraft
29, 486–488. doi: 10.1115/1.3640593 drag modelling. Progr. Aerosp. Sci. 64, 56–65. doi: 10.1016/j.paerosci.2013.
Guettler, D. B. (2007). Satellite Attitude Control Using Atmospheric Drag. Master’s 09.001
thesis, Air Force Institute of Technology. Psiaki, M. L. (2004). Nanosatellite attitude stabilization using passive
Halsmer, D. M., and Mingori, D. L. (1995). Nutational stability and passive control aerodynamics and active magnetic torquing. J. Guid. Control Dyn. 27,
of spinning rockets with internal mass motion. J. Guid. Control Dyn. 18, 347–355. doi: 10.2514/1.1993
1197–1203. doi: 10.2514/3.21525 Scholz, C., Romagnoli, D., Dachwald, B., and Theil, S. (2011). Performance analysis
Hamidi-Hashemi, H. (1993). “Liapunov analysis of a two dimensional of an attitude control system for solar sails using sliding masses. Adv. Space Res.
unconstrained particle motion in a rigid body spinning about the thrust axis,” 48, 1822–1835. doi: 10.1016/j.asr.2011.05.032
in Circuits and Systems, 1993., Proceedings of the 36th Midwest Symposium on Sentman, L. (1961). Free Molecule Flow Theory and Its Application to the
(Detroit, MI), Vol.2. 971–973. Determination of Aerodynamic Forces. Sunnyvale, CA: Lockheed Missile and
Hedin, A. E., Fleming, E. L., Manson, A. H., Schmidlin, F. J., Avery, S. K., Clark, Space Co.
R. R., et al. (1996). Empirical wind model for the middle and lower atmosphere. Shao, A., Koltz, E. A., and Wertz, J. R. (2014). “Quantifying the cost reduction
J. Atmos. Terrestr. Phys. 58, 1421–1447. doi: 10.1016/0021-9169(95)00122-0 potential for earth observation satellites,” in 12th Reinventing Space Conference,
ISO 14222 (2013). ISO 14222 Space Environment (Natural and Artificial). Earth 18-20 November 2014 (London).
Upper Atmosphere. Technical Report ISO 14222:2013, ISO. Sutton, E. K. (2009). Normalized force coefficients for satellites with elongated
Janssens, F. L., and van der Ha, J. C. (2014). Stability of spinning satellite under shapes. J. Spacecr. Rockets 46, 112–116. doi: 10.2514/1.40940
axial thrust, internal mass motion, and damping. J. Guid. Control Dyn. 38, Virgili, J., Roberts, P. C. E., and Hara, N. C. (2015). Atmospheric interface reentry
761–771. doi: 10.2514/1.G000123 point targeting using aerodynamic drag control. J. Guid. Control Dyn. 38,
Killeen, T. L., Hays, P. B., Spencer, N. W., and Wharton, L. E. (1982). Neutral winds 403–413. doi: 10.2514/1.G000884
in the polar thermosphere as measured from dynamics explorer. Geophys. Res. Virgili-Llop, J. (2014). Spacecraft Flight in the Atmosphere. Ph.D. thesis, School of
Lett. 9, 957–960. doi: 10.1029/GL009i009p00957 Engineering, Cranfield University.
King-Hele, D. (1992). The upper atmosphere as sensed by satellite orbits. Planet. Virgili-Llop, J., Polat, H. C., and Romano, M. (2016). Using shifting masses to reject
Space Sci. 40, 223–233. doi: 10.1016/0032-0633(92)90060-2 aerodynamic perturbations and to maintain a stable attitude in very low earth
King-Hele, D., and Walker, D. M. (1988). Upper-atmosphere zonal winds orbit. Adv. Astron. Sci. 158, 2129–2148.
from satellite orbit analysis: an update. Planet. Space Sci. 36, 1085–1093. Virgili-Llop, J., Roberts, P. C., Hao, Z., Tomas, L. R., and Beauplet, V. (2014a).
doi: 10.1016/0032-0633(88)90062-1 “Very low earth orbit mission concepts for earth observation. benefits and
King-Hele, D. G. (1987). Satellite Orbits in an Atmosphere: Theory And Application. challenges,” in Reinventing Space Conference, 18-21 November, Vol. BIS-RS-
London: Blackie Academic and Professional. 2014-37 (London).

Frontiers in Robotics and AI | www.frontiersin.org 18 February 2019 | Volume 6 | Article 7

Virgili-Llop et al. Attitude Stabilization by CoM Shifting

Virgili-Llop, J., Roberts, P. C. E., Palmer, K., Hobbs, S., and Kingston, J. (2014b). Wie, B., Weiss, H., and Arapostathis, A. (1989). Quarternion feedback regulator
Descending sun-synchronous orbits with aerodynamic inclination correction. for spacecraft eigenaxis rotations. J. Guid. Control Dyn. 12, 375–380.
J. Guid. Control Dyn. 38, 831–842. doi: 10.2514/1.G000183
Wertz, J., Sarzi-Amade, N., Shao, A., Taylor, C., and Van Allen, R. (2012). Conflict of Interest Statement: The authors declare that the research was
“Moderately elliptical very low orbits (mevlos) as a long-term solution to conducted in the absence of any commercial or financial relationships that could
orbital debris,” in 26th Annual AIAA/USU Conference on Small Satellites, Vol. be construed as a potential conflict of interest.
SSC12-IV-6 (Logan, UT).
Wie, B. (2004). Solar sail attitude control and dynamics, part 1. J. Guid. Control Copyright © 2019 Virgili-Llop, Polat and Romano. This is an open-access article
Dyn. 27, 526–535. doi: 10.2514/1.11134 distributed under the terms of the Creative Commons Attribution License (CC BY).
Wie, B., and Barba, P. M. (1985). Quaternion feedback for spacecraft large angle The use, distribution or reproduction in other forums is permitted, provided the
maneuvers. J. Guid. Control Dyn. 8, 360–365. doi: 10.2514/3.19988 original author(s) and the copyright owner(s) are credited and that the original
Wie, B., and Murphy, D. (2007). Solar-sail attitude control design for a publication in this journal is cited, in accordance with accepted academic practice.
flight validation mission. J. Spacecr. Rockets 44, 809–821. doi: 10.2514/1. No use, distribution or reproduction is permitted which does not comply with these
22996 terms.

Frontiers in Robotics and AI | www.frontiersin.org 19 February 2019 | Volume 6 | Article 7