Algorithms For Small Satellite Formation Flying

Air Force Institute of Technology
AFIT Scholar
Theses and Dissertations Student Graduate Works
3-22-2018
Algorithms for Small Satellite Formation Flying

Robert B. LaRue
Follow this and additional works at: https://scholar.afit.edu/etd
Part of the Astrodynamics Commons, and the Navigation, Guidance, Control and Dynamics Commons
Recommended Citation
LaRue, Robert B., "Algorithms for Small Satellite Formation Flying" (2018). Theses and Dissertations.
1776.
https://scholar.afit.edu/etd/1776
This Thesis is brought to you for free and open access by the Student Graduate Works at AFIT Scholar. It has been
accepted for inclusion in Theses and Dissertations by an authorized administrator of AFIT Scholar. For more
information, please contact AFIT.ENWL.Repository@us.af.mil.
Algorithms for Small
Satellite Formation Flying
THESIS
Robert B. LaRue, 2nd Lieutenant, USAF
AFIT-ENY-MS-18-M-273
DEPARTMENT OF THE AIR FORCE

AIR UNIVERSITY
AIR FORCE INSTITUTE OF TECHNOLOGY

Wright-Patterson Air Force Base, Ohio
DISTRIBUTION A. APPROVED FOR PUBLIC RELEASE; DISTRIBUTION IS UNLIMITED.

The views expressed in this thesis are those of the author and do not reflect the
official policy or position of the United States Air Force, Department of Defense, or
the United States Government.

THESIS
Presented to the Faculty

Department of Aeronautics and Astronautics
Graduate School of Engineering and Management
Air University
Air Education and Training Command
In Partial Fulfillment of the Requirements for the
Degree of Master of Science in Astronautical Engineering
Robert B. LaRue, B.S.

2nd Lieutenant, USAF
March 2018
DISTRIBUTION A. APPROVED FOR PUBLIC RELEASE; DISTRIBUTION IS UNLIMITED.

ALGORITHMS FOR SMALL

SATELLITE FORMATION FLYING
Robert B. LaRue, B.S.

2nd Lieutenant, USAF
Committee Membership:
Lt Col Kirk W. Johnson, PhD

Chair
William E. Wiesel, PhD

Member
Capt Joshuah A. Hess, PhD

Member
Abstract
This thesis presents algorithms for spacecraft formation flying using impulsive-
thrust and low-thrust methods. The general circular orbit formation initial conditions
are derived in terms of equinoctial elements. Physical significance of the bounded rel-
ative motion parameters is presented for the case of general circular orbits. The
developed algorithms are posed in terms of equinoctial elements for a singularity-free
approach. The algorithms are assessed by numerical propagation of the inertial equa-
tions of motion with J2 and drag perturbations. Methods are presented for minimizing
the ∆V required for formation initialization. An examination of the performance of
open-loop and closed-loop control is provided for formation initialization and recon-
figuration. The effects of differential drag on small satellite formations is analyzed.
The developed algorithms are used to examine the trade space and quantify how
spacecraft design parameters affect formation flying scenarios.
iv
Acknowledgements
I would like to thank my advisor, Lieutenant Colonel Kirk Johnson, for his mentorship
and guidance throughout my time here at AFIT. Completing my research was made
possible by drawing on his expertise and knowledge. Additional thanks go to the
student and faculty members of the RPO working group, whose weekly feedback
greatly improved the quality and clarity of my work.
I would also like to thank Dr. William Wiesel and Captain Joshuah Hess for
serving on my thesis committee. I’m also very grateful to the Air Force Research
Laboratory for supporting my research here at AFIT. Finally, I would like to thank
the U.S. Air Force for allowing me to spend my first assignment as a full-time graduate
student. It has been an amazing opportunity and something for which I will be
eternally grateful.
Robert B. LaRue
v
Table of Contents
Page
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
I. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
II. Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1 Relative Motion and Formation Flying . . . . . . . . . . 4
2.2 Satellite Formations . . . . . . . . . . . . . . . . . . . . 6
2.2.1 PCO and GCO Formations . . . . . . . . . . . 6
2.2.2 Along-Track Orbits . . . . . . . . . . . . . . . . 9
2.3 Orbit Element Sets . . . . . . . . . . . . . . . . . . . . . 10
2.4 Orbit Element Differences . . . . . . . . . . . . . . . . . 13
2.5 Orbital Perturbations . . . . . . . . . . . . . . . . . . . 14
2.5.1 Perturbations due to an Aspherical Gravity Field 15
2.5.2 Perturbations due to Atmospheric Drag . . . . . 16
2.6 Mean Orbit Elements . . . . . . . . . . . . . . . . . . . 18
2.7 Small Satellites . . . . . . . . . . . . . . . . . . . . . . . 21
2.8 State Transition Matrix . . . . . . . . . . . . . . . . . . 22
2.9 Maneuvering Techniques . . . . . . . . . . . . . . . . . . 23
2.9.1 Impulsive Maneuvering . . . . . . . . . . . . . . 24
2.9.2 Continuous-Thrust Maneuvering . . . . . . . . . 25
2.10 Propulsion Technologies . . . . . . . . . . . . . . . . . . 26
2.10.1 Chemical Propulsion . . . . . . . . . . . . . . . 27
2.10.2 Electrical Propulsion . . . . . . . . . . . . . . . 27
III. Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.1 Formation Design . . . . . . . . . . . . . . . . . . . . . . 29
3.2 Impulsive Formation Control . . . . . . . . . . . . . . . 34
3.2.1 Impulsive Maneuvering Algorithm . . . . . . . . 35
3.2.2 Impulsive Formation Maintenance Algorithm . . 37
3.3 Continuous-Thrust Formation Control . . . . . . . . . . 39
3.4 Coordinate Transformations . . . . . . . . . . . . . . . . 41
3.4.1 ECI to Hill Frame Transformation . . . . . . . . 41
vi
Page
3.4.2 Hill Frame to ECI Transformation . . . . . . . . 42
3.4.3 Determining ECI Position and Velocity from Equinoc-
tial Elements . . . . . . . . . . . . . . . . . . . 45
3.5 Initial Conditions for Formation Initialization . . . . . . 47
3.6 Numerical Propagation of Solution . . . . . . . . . . . . 57
3.7 Algorithm Overview . . . . . . . . . . . . . . . . . . . . 59
IV. Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.1 Formation Initialization . . . . . . . . . . . . . . . . . . 61
4.1.1 Impulsive-Thrust Initialization . . . . . . . . . . 61
4.1.2 Low-Thrust Initialization . . . . . . . . . . . . . 63
4.2 Formation Reconfiguration . . . . . . . . . . . . . . . . . 68
4.2.1 Impulsive-Thrust Reconfiguration . . . . . . . . 68
4.2.2 Low-Thrust Reconfiguration . . . . . . . . . . . 75
4.3 Drag Considerations . . . . . . . . . . . . . . . . . . . . 79
V. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.1 Contributions and Key Findings . . . . . . . . . . . . . . 84
5.2 Limitations of the Algorithms . . . . . . . . . . . . . . . 85
5.3 Recommendations for Future Work . . . . . . . . . . . . 86
Appendix A. State Transition Matrix . . . . . . . . . . . . . . . . . . 88
Appendix B. Additional Results . . . . . . . . . . . . . . . . . . . . . 90
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
vii
List of Figures
Figure Page
1. Hill Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2. PCO Relative Orbit Geometry at t = 0 . . . . . . . . . . . . . 8
3. GCO Relative Orbit Geometry at t = 0 . . . . . . . . . . . . . 9
4. Orbit Orientation . . . . . . . . . . . . . . . . . . . . . . . . . 11
5. Categories of Small Satellites . . . . . . . . . . . . . . . . . . . 21
6. Impulsive Maneuver . . . . . . . . . . . . . . . . . . . . . . . . 25
7. Spacecraft Velocity Relative to ECI Frame . . . . . . . . . . . 43
8. CSD Ejection Velocity . . . . . . . . . . . . . . . . . . . . . . . 47
9. PCO Initialization Maneuver With N = 2, ∆tmax = 3 Orbits . 62
10. PCO Initialization Maneuver With N = 2, ∆tmax = 1 Orbit . . 63
14. PCO Initialization Maneuver With N = 5, ∆tmax = 8 Orbits,
Closed-Loop Feedback . . . . . . . . . . . . . . . . . . . . . . . 65
15. Low-Thrust PCO Initialization Maneuver, ∆t = 4.5 Orbits . . 66
16. Low-Thrust PCO Initialization Maneuver, ∆t = 8 Orbits . . . 67
17. Low-Thrust PCO Initialization Maneuver, ∆t = 20 Orbits . . . 67
18. 1 km → 2 km GCO Reconfiguration with ∆tmax = 3 orbits, N = 2 69
22. 1 km → 2 km GCO Reconfiguration with ∆tmax = 5 orbits,
N =5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
viii
Figure Page
25. GCO Formation Reconfiguration Results Using Impulsive-Thrust 72
26. Approximate Linear Relationship For Impulsive GCO Reconfig-
uration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
27. 1 km → 2 km ATO Reconfiguration with ∆tmax = 3 orbits, N = 2 74
29. 1 km → 2 km ATO Reconfiguration r with ∆tmax = 3 orbits,
N =4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
32. 1 km → 2 km Low-thrust GCO Reconfiguration Maneuver . . . 76
33. 1 km → 3 km Low-Thrust GCO Reconfiguration Maneuver . . 77
34. 1 km → 2 km Low-Thrust GCO Reconfiguration Maneuver with
500 kg, 23 mN Thrust Satellite . . . . . . . . . . . . . . . . . . 77
10 kg, 4 mN Thrust Satellite . . . . . . . . . . . . . . . . . . . 78
37. 0.5 km → 1 km Low-Thrust GCO Reconfiguration Maneuver . 79
38. 1 km PCO Initialization Maneuver with Worst-Case Differential
Drag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
Drag, 6 Impulses . . . . . . . . . . . . . . . . . . . . . . . . . . 81
Drag using Low-Thrust, ∆t = 4.5 orbits . . . . . . . . . . . . . 82
Drag using Low-Thrust, ∆t = 8 orbits . . . . . . . . . . . . . . 82
Drag using Low-Thrust, ∆t = 20 orbits . . . . . . . . . . . . . 83
43. 1 km → 2 km PCO Reconfiguration with ∆tmax = 3 orbits, N = 2 90
ix
Figure Page
46. 1 km → 2 km PCO Reconfiguration with ∆tmax = 3 orbits,
N =5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
48. 1 km → 2 km Low-Thrust PCO Reconfiguration Maneuver . . 92
49. 1 km → 3 km Low-Thrust PCO Reconfiguration Maneuver . . 93
50. 1 km → 2 km Low-Thrust PCO Reconfiguration Maneuver with
500 kg, 23 mN Thrust Satellite . . . . . . . . . . . . . . . . . . 93
53. 0.5 km → 1 km Low-Thrust PCO Reconfiguration Maneuver . 95
x
List of Tables
Table Page
1. Exponential Density Model . . . . . . . . . . . . . . . . . . . . 17

2. Summary of PCO and GCO Initial Conditions . . . . . . . . . 34
3. Spacecraft Parameters . . . . . . . . . . . . . . . . . . . . . . . 66
4. 1 km → 2 km PCO Reconfiguration Maneuver Results . . . . . 69
5. PCO Reconfiguration Using Low Thrust . . . . . . . . . . . . . 79
xi
I. Introduction
With the miniaturization of computers and satellite components, the achievable
performance of small satellites (SmallSats) is rapidly increasing. Such technological
developments make it possible for SmallSats to perform missions that previously re-
quired larger and more expensive satellites. Greater access to space at lower cost will
make it possible for universities and small companies to launch satellite networks and
space exploration missions. Cube satellites (CubeSats) in particular can be devel-
oped rapidly for low cost, and have seen ever-increasing use since their introduction
by Stanford and Cal-Poly in 1999 [1] [2]. Originally proposed as a low-cost method
to get students involved in space, CubeSats have evolved into a flourishing area of
research that has potential for a variety of missions. The low-cost access to space
that CubeSats provide has intrigued major space companies and government admin-
istrations as a way to demonstrate new technologies or perform missions that are not
feasible with larger satellites [3].
A satellite formation is defined by NASA’s Goddard Space Flight Center (GSFC)

as “The tracking or maintenance of a desired relative separation, orientation or posi-
tion between or among spacecraft” [4]. Formation flying has many potential applica-
tions in which a formation of spacecraft work together to accomplish a task that would
otherwise require a significantly larger or more complex spacecraft. One such appli-
cation is distributed sensing, where a satellite formation could achieve performance
that would conventionally require a massive radar or telescope. Satellite formations
also allow for greater reliability and resiliency, as it would be much easier to replace a
malfunctioning satellite than it would be to repair a component on a larger satellite.
1
A number of formation flying missions have been conducted to date, such as
the Magnetospheric Multiscale (MMS) mission [5] [6], and the Prototype Research
Experiments and Space Mission Technology Advancement (PRISMA) [7], among oth-
ers. There have also been CubeSat formation flying missions, such as AeroCube-4 [8],
which demonstrated formation flying by changing the satellites’ drag coefficient to
reconfigure the formation over a period of weeks. CubeSat Proximity Operations
Demonstration (CPOD) [9] aims to demonstrate proximity operations and docking
using two 3U CubeSats. The CanX-4 and CanX-5 spacecraft demonstrated precision
formation flying using picosatellites, including projected circular orbit (PCO) forma-
tions [10]. An overview of many current and future CubeSat formation flying missions
can be found in [8].
While there are benefits to using a formation of satellites as opposed to a single

spacecraft, there are several technical difficulties that arise. Accurate control of the
formation geometry is a necessity for mission success in many cases. For applications
such as distributed sensing, where the spacecraft in the formation work together to
perform the function of a radar or telescope, it may be desirable to resize the formation
periodically to change the resolution of the distributed system. This is an example of
a formation reconfiguration problem.
There have been numerous studies on the reconfiguration of satellite relative

orbits. Vaddi et al. developed an analytical two-impulse control scheme to transfer
a deputy between given initial and final configurations [11], including unperturbed
PCO and GCO (General Circular Orbit) formations. Palmer investigated relative
low-thrust transfers for Keplerian orbits [12]. Vignal and Pernicka examined the per-
formance of linear vs. nonlinear controllers for low-thrust formation keeping [13]. Yan
and Alfriend used the Gim-Alfriend state transition matrix [14] to analytically derive
low-thrust control laws for perturbed relative motion [15]. Cho and Park developed
a computationally efficient solution to the reconfiguration problem that does not re-
quire inverting the fundamental dynamics matrix [16]. Lee and Park provided an
analytical solution to the low-thrust formation reconfiguration problem including dif-
2
ferential gravity, J2 , and chief eccentricity [17]. Massari and Bernelli-Zazzera derived
optimal low-thrust control laws subject to multiple constraints using the interior-
point method [18]. Cho, Park, Yoo, and Choi derived an analytical solution for
formation reconfigurations on elliptical reference orbits using the Tschauner-Hempel
equations [19] [20]. These and other studies have provided crucial insight into the
dynamics and control of the formation reconfiguration problem.
The work presented in this thesis will leverage the literature to build algorithms
for formation initialization, reconfiguration, and maintenance. These algorithms will
be used to analyze a variety of scenarios involving small satellite formation flying.
The goal is to examine the tradespace and determine how spacecraft characteristics
such as thrust, mass, and specific impulse affect formation flying maneuvers. Initial
conditions will be chosen to simulate realistic scenarios [21] [22]. The effectiveness
of the guidance and control algorithms will be assessed by a numerical propagation
using absolute equations of motion with J2 and drag effects included, to capture the
full nonlinear equations of motion and major perturbations.
To summarize, the goals of the research presented are as follows:
1. Develop an algorithm for formation flying that accounts for situations unique
to SmallSats
2. Compare and evaluate a variety of formation types
3. Examine design tradespace to see how spacecraft parameters affect formation

flying scenarios
This Introduction chapter will be followed by a Background chapter, which sum-

marizes the concepts that form the foundation for the work presented in this thesis.
Next is the Methodology chapter, which will outline the developed algorithms and
present the overall process for obtaining solutions. The Results chapter shows simula-
tion outcomes to quantify the performance of the algorithms. Finally, the Conclusion
will summarize the Results and state the significant outcomes of the research.
3
II. Background
This chapter will review the basic principles that form the foundation for the research
presented in this thesis.
2.1 Relative Motion and Formation Flying
One of original models for satellite relative motion is the Hill/Clohessy-Wiltshire

(HCW) equations [23] [24]. These equations were originally derived for rendezvous
and docking in the early space program. The homogenous form of the HCW equations
is provided below:
ẍ − 2nẏ − 3n2 x = 0 (1)
ÿ + 2nẋ = 0 (2)
z̈ + n2 z = 0 (3)
where n is the chief satellite mean motion and x, y, and z are the Hill frame compo-
nents (Fig. 1). These are linearized equations of motion that make several assump-
tions, such as zero chief eccentricity, no perturbation forces, and that the deputy
satellite is close to the chief (d << r, where d is the distance from the chief to the
deputy, and r is the orbit radius of the chief). Despite the assumptions, the HCW
equations are an attractive option because they are a set of autonomous equations
of motion (they do not explicitly depend on time) with a simple closed form solu-
tion [4] [25]:
2ẏ0 ẋ0 2ẏ0

x(t) = 4x0 + + sin (nt) − 3x0 + cos (nt) (4)
n n n
2ẋ0 4ẏ0 2ẋ0

y(t) = −(6nx0 + 3ẏ0 )t + cos (nt) + 6x0 + sin (nt) − + y0 (5)
n n n
ż0
z(t) = z0 cos (nt) + sin (nt) (6)
n
where the subscript 0 refers to the initial conditions.
4
Figure 1: Hill Frame
The HCW equations are typically an acceptable model for modeling relative mo-
tion in cases where the chief orbit is very close to circular, the deputy is within a few
kilometers of the chief, and forces other than gravity remain small. The inaccuracy of
the HCW model will rapidly increase with the degree to which these assumptions are
violated. Even in cases where these assumptions hold, predicting motion over long
periods of time using the HCW model will be inaccurate due to neglecting pertur-
bation forces, particularly J2 . There have been numerous efforts in the literature to
develop higher fidelity relative motion models that relax the simplifying assumptions
used in deriving the HCW equations. For example, the Schweighart-Sedwick model is
a linearized relative motion model similar to HCW, but includes J2 perturbations [26].
Gim and Alfriend developed a state transition matrix (STM) that accommodates an
elliptical chief orbit and J2 perturbations using the geometric method, and provides
transformations between Hill components and orbit element differences [14]. There
are also models such as the Tschauner-Hempel equations [20] that accommodate chief
eccentricity by using true anomaly rather than time as the independent variable. For
a thorough overview and comparison of many of the available relative motion models,
the reader is referred to [27].
5
Even though the HCW equations leverage a variety of assumptions, they can
still be used to reveal fundamental characteristics about satellite relative motion.
Equation (6) demonstrates that the cross-track motion is a simple harmonic. This
cross-track oscillation can be negated by choosing z0 = ż0 = 0. The first term in Eq.
(5) is often referred to as the secular term, as it is the only term in the HCW solution
that grows with time (all other t terms are periodic). The secular term results in
general instability for in-plane (x-y plane) motion. However, a stable subspace exists
under certain conditions where the secular term is eliminated. The initial conditions
that eliminate the secular drift are:
ẏ0 = −2nx0 (7)
Many formation flying strategies include choosing initial conditions such that Eq.
(7) is satisfied, therefore eliminating drift and bounding the motion of the deputy to
the chief. It is important to note that initial conditions that satisfy Eq. (7) bound
motion only for the HCW model (general bounded motion occurs when the spacecraft
in a formation have equal energies). If a higher fidelity model is being used, there
are additional constraints that must be used to fully eliminate drift [4]. Stability
in the context of the HCW model is local – there are initial conditions that violate
Eq. (7) but satisfy global boundedness. In spite of these caveats, there are interesting
geometric properties that emerge from the HCW solution. These geometric properties
form the basis for PCO and GCO orbits.
2.2 Satellite Formations
2.2.1 PCO and GCO Formations. This section will describe the basic
characteristics of two particular types of satellite formation - projected circular orbits
(PCOs) and general circular orbits (GCOs). A PCO is a trajectory that forms a
circular projection onto the y-z (horizontal) plane of the Hill frame. These orbits
have a variety of applications for Earth-observing missions, as this formation leads to
6
equal spacing of the deputy spacecraft in the horizontal plane. A GCO is a trajectory
that forms a three-dimensional circle in the Hill frame. This formation type may
be suitable for applications where constant spacing between the chief and deputy
satellites is desired. Both formations are defined in terms of the phase-magnitude
form of the stable subspace HCW solution [4]:
x(t) = ρx sin (nt + αx ) (8)
y(t) = ρy + 2ρx cos (nt + αx ) (9)
z(t) = ρz sin (nt + αz ) (10)
These equations are a form of the HCW solution with the initial conditions chosen
such that Eq. (7) is satisfied. The variables ρx , ρy , and ρz are magnitude (distance)
parameters, and the variables αx and αz are phase (angle) parameters. Equations (8)
and (9) show a fundamental property of bounded relative motion – the 2:1 in-plane
ellipse. The expressions for the bounded relative motion parameters in terms of Hill
frame components are presented below:
p
ẋ20 + x20 n2
ρx = (11)
n
ẋ0
ρy = y0 − 2 (12)
n
p
ż02 + z02 n2
ρz = (13)
n
αx = atan2(nx0 , ẋ0 ) (14)
αz = atan2(nz0 , ż0 ) (15)
where atan2 is the quadrant-specific inverse tangent function. In the special case
√
where αx = αz and ρz = 3ρx , the relative orbit is a three-dimensional circle with
radius 2ρx (a GCO) centered at (0, ρy , 0). In the case where αx = αz and ρz = 2ρx ,
the relative orbit is a PCO centered at (0, ρy , 0). The PCO and GCO constraints
7
are typically defined with ρy = 0, to center the formation on the chief. The angle
parameters generally do not have any immediate physical significance. However, in
the special case of a PCO, αz is the phase of the deputy satellite in the projected
circle (measured from the +y axis) at the time of the chief satellite’s equator crossing.
Physical significance of the bounded relative motion parameters for a PCO are shown
in Fig. 2 [28]. Physical significance of the bounded relative motion parameters for a
GCO are shown in Fig. 3.
Figure 2: PCO Relative Orbit Geometry at t = 0
There are other special cases of bounded relative motion in the literature, but
PCOs and GCOs have been the subject of much of the current work. PCOs have been
8
Figure 3: GCO Relative Orbit Geometry at t = 0
a particular area of interest, due to their attractive properties and J2 -invariant nature.
Natural perturbations cause a phase rotation α̇ in the PCO formation that allow it
to be controlled for long periods of time with little fuel cost , as well as balancing fuel
consumption of the various satellites in the formation [29] [30]. Much of the analysis
presented in the results section will be focused on PCO and GCO formations due to
their unique properties and variety of applications.
2.2.2 Along-Track Orbits. Along-track orbits (ATOs) are a type of space-

craft formation where the deputy satellite has an offset in the y-direction of the Hill
frame. ATOs are sometimes referred to as “leader-follower” formations. The initial
9
conditions for an along-track orbit can be characterized as:
   
x 0
 0  
 y0  yd 
   
   
   
 z0   0 
X̄0 = 
 = 
   (16)
ẋ0   0 
   
   
 ẏ0   0 
   
ż0 0
where yd is the desired offset in the y-direction. In the HCW model, an ATO will
result in a relative trajectory with y(t) = y0 , and x(t) = z(t) = 0 (observe Eqs. (4 -
6)). As a result, the relative trajectory will appear to be stationary in the Hill frame.
In reality, orbital perturbations and nonlinearities will cause the trajectory to deviate
somewhat from the (0, y0 , 0) point.
There are other variants of the ATO where one or more of the deputy satellites
are given an offset in the z-direction. In this case, the satellite will have oscillatory
cross-track motion, but will still maintain a constant offset in the y-direction (since
cross-track motion is decoupled from in-plane motion for the HCW model). This
type of modified ATO may be desirable in cases when cross-track separation between
satellites is needed. It is worth noting that since the z-motion is oscillatory, the
desired cross-track separation would only occur at one point in the orbit.
2.3 Orbit Element Sets
There are a variety of available orbit element sets that may be used to charac-
terize a satellite orbit. The classical orbit elements are defined as:
¯ = (a, e, i, h, g, l)
cl (17)
where a is the semimajor axis, e is the eccentricity, i is the inclination, h is the

right ascension of the ascending node, g is the argument of perigee, and l is the mean
10
Figure 4: Orbit Orientation
anomaly. The orbit elements a and e define the size and shape of the orbit, while the
elements i, h, and g define the orientation of the orbit in space. The mean anomaly
l defines the position of the satellite in the orbit. Together, these six elements fully
define the trajectory of a satellite. Figure 4 shows the orientation of an orbit in space.
X̂, Ŷ , and Ẑ are respectively the 1-, 2-, and 3-direction of the Earth-centered inertial
(ECI) frame. X̂ is known as the vernal equinox, or the first point of Aries. Ẑ points
through the north pole, and Ŷ = Ẑ × X̂.
The right ascension of the ascending node (h), sometimes simply referred to
as the right ascension or the node, is the angle between X̂ and the nodal vector ~n.
The nodal vector locates the point where the satellite orbit crosses the X̂-Ŷ plane
from the −Ẑ to +Ẑ direction. g measures the angle from ~n to the orbit periapsis.
i measures the tilt of the orbit plane, and is measured between K̂ and the orbit
angular momentum vector. The inclination can lie between 0◦ and 180◦ . Orbits with
inclinations between 0◦ and 90◦ are referred to as prograde orbits, as the satellite
moves through its orbit in the same direction as the Earth’s rotation. Alternatively,
11
orbits with inclinations between 90◦ and 180◦ are retrograde orbits, as the satellite
moves through its orbit in the opposite direction of the Earth’s rotation. Orbits with
inclinations of 0◦ or 180◦ are equatorial orbits, the former being prograde and the
latter retrograde. An inclination of 90◦ corresponds to a polar orbit.
The true anomaly ν measures the angle between the satellite and the orbit
periapsis. The mean anomaly l is a fictitious angle that measures the angle to the
satellite if it moved in a circular orbit with an equal period to the real orbit, and is
related to the true anomaly by Kepler’s equation. While l may be less intuitive than
ν, it has the advantage of a constant rate (the mean motion) in the two-body problem
for any elliptical orbit. Either l or ν may be used to define the location of a satellite
in its orbit. For the special case of a circular orbit, l and ν are equivalent. For the
interested reader who wishes to learn more about the classical orbit elements and
Kepler’s equation, there are thorough discussions of these topics in refs. [25] and [31].
The classical orbit elements have the advantage of using parameters that all have
intuitive significance. However, they are singular for orbits that are circular (e = 0) or
equatorial (i = 0), as g is undefined for the case of circular orbits and h is undefined for
equatorial orbits. These singularities make using classical orbit elements impractical,
as satellites frequently are placed in orbits where they arise. Geostationary spacecraft
are a common example of satellites with orbits that are both circular and equatorial.
The nonsingular elements are an alternative orbit element set that can accomo-
date circular orbits. These elements are defined as:
ns
¯ = (a, λ, i, q1 , q2 , h) (18)
where λ = l + g, q1 = e cos g, and q2 = e sin g. It is important to note that the

nonsingular elements are still singular for equatorial orbits, as h is undefined in this
case. In order to accommodate orbits that are both circular and equatorial, the
12
equinoctial elements may be used. The equinoctial elements are defined as:
eq
¯ = (a, Λ, q̃1 , q̃2 , p1 , p2 ) (19)
where Λ = l + g + h, q̃1 = e cos (g + h), q̃2 = e sin (g + h), p1 = tan (i/2) cos h,
and p2 = tan (i/2) sin h. The equinoctial elements have a singularity for retrograde
equatorial orbits (i = 180 deg), but there are alternative formulations that account for
this case [32]. However, an Earth satellite has not yet been launched into this orbit
[33], so the equinoctial elements are typically considered a singularity-free element
set. The notation for the equinoctial elements is not universally standardized, so it
is important to understand which definition is being used for each case. Equinoctial
elements have the advantage of the least amount of singularities of the presented
element sets, but the disadvantage of orbit elements whose physical significance is
not intuitive. In general, the equinoctial elements are more difficult to work with
than classical or nonsingular elements, as they tend to lead to complicated algebraic
expressions. It is for this reason that nonsingular elements are generally preferable to
equinoctial elements for the case of inclined circular orbits.
Three sets of orbit elements have been presented here, but there are a multitude
of orbit element sets available in the literature, each with their own unique set of
advantages and disadvantages. For a review of the available orbit element sets, the
reader is referred to refs. [34] and [28].
2.4 Orbit Element Differences
The relative orbit of a satellite can be characterized using orbit element differ-
ences:
δ oe ¯ d − oe
¯ = oe ¯c (20)
where oe
¯ d are the orbit elements of the deputy and oe
¯ c are the orbit elements of
the chief. Note that oe ¯ ns,
¯ can be any of the presented orbit element sets (cl, ¯ or
eq).
¯ Given the chief orbit elements oe
¯ c and the element differences, the position of
13
the deputy at any time can be solved for using Kepler’s equation (or an appropriate
modified version). Equation (20) makes no assumptions about the chief satellite orbit,
nor does it require that the deputy be in close proximity to the chief. The Hill frame
state vector has six fast variables (variables that change with time), while the orbit
element sets only have one (l, λ, or Λ) in the two-body problem, thus simplifying the
the satellite position computation.
There are inferences that can be made about a relative trajectory described
using orbit element differences that require no computation. The δi and δh of a
deputy satellite determine the magnitude of its cross-track motion. Differences in
semimajor axis (δa) result in the chief and deputy satellite having a different period,
causing them to drift apart [35]. For the case of the two-body problem, a δa of zero
will result in no relative drift. It is for this reason that the equation of constraint
δa = 0 is equivalent to Eq. (7). However, a modified version of this constraint must
be chosen in the presence of perturbations.
2.5 Orbital Perturbations
The motion of a spacecraft is primarily governed by the gravity of the primary

body which it orbits. This motion can be described using the two-body equation:
µ
~r¨ = − 3 ~r (21)
r
where µ is the gravity parameter of the primary body, and ~r is the vector from the
center of the primary body to the spacecraft (r = ||~r||). Equation (21) relies on
several assumptions, namely that the mass of the spacecraft is negligible compared
to the central body, the chosen coordinate system is inertial, the primary body is
spherically symmetric with uniform density, and no forces act on the system other
than gravity between the spacecraft and the primary body. For the case of Earth-
orbiting spacecraft, the first two of these assumptions can generally be assumed to be
true. However, the latter two are not – the Earth is not a perfect sphere, and there
14
are a variety of forces that act on the system other than the force of gravity between
the two bodies. As long as these forces are small compared to the two-body gravity,
they can be modeled as perturbations (small deviations) from two-body gravity. In
this case, equation 21 can be modified:
µ
~r¨ = − 3 ~r + ~ad (22)
r
where ~ad are the disturbances (perturbations) from the two-body trajectory. Exam-
ples of typical perturbations are an aspherical gravity field (J2 ), atmospheric drag,
third-body gravity (sometimes called lunisolar perturbations in Earth orbit), and
solar radiation pressure. These perturbations are present for all Earth-orbiting space-
craft, but the importance of each is varied depending the spacecraft’s orbit regime.
In low-Earth orbit (LEO), the largest perturbations are due to aspherical gravity
and atmospheric drag. At higher altitudes such as geosynchronous orbits (GEO),
atmospheric drag is not a major concern, but solar radiation pressure and third-body
gravity from the Sun and Moon have a significant effect.
2.5.1 Perturbations due to an Aspherical Gravity Field. The analysis pre-

sented in this thesis is focused on SmallSats in low-Earth orbit (LEO). In LEO, J2
is typically the largest perturbation, followed by atmospheric drag. In many cases,
modeling only J2 perturbations is a reasonable approximation of reality for spacecraft
in LEO. J2 is a zonal harmonic term due to the Earth’s equatorial bulge, and is the
largest aspherical gravity perturbation term (by several orders of magnitude). The
acceleration component due to J2 can be modeled as [36]:
 X 
1 − 5( Zr )2
r
3 µ Re 2  
~aJ2 = − J2 2 Z 2 Y 
1 − 5( r ) r  (23)

2 r r

 
Z 2 Z
3 − 5( r ) r
where Re is the radius of the Earth and X, Y , and Z are the position components
in the Earth-centered inertial (ECI) reference frame (Fig. 4). Equation (23) can be
15
inserted as part of the ~ad term in Eq. (22) to account for J2 perturbations. In certain
cases where a high degree of accuracy is desired, it may be necessary to augment Eq.
(23) to account for additional spherical harmonics.
2.5.2 Perturbations due to Atmospheric Drag. The perturbation due to

drag can be modeled as:
1 CD S
~aaero = − ζ ||~vrel ||~vrel (24)
2 m
where ζ is the density of the atmosphere, CD is the spacecraft drag coefficient, S is
the projected area (surface area of spacecraft exposed to drag), m is the spacecraft
mass, and ~vrel is the velocity relative to the atmosphere. CD is determined by the
spacecraft design, but can typically be approximated as CD = 2.2 [25]. S is a function
of both the spacecraft design and its attitude (how it is oriented in space). ~vrel can
be found in terms of ECI components using the following equation:
~vrel = ~v − (~ωE × ~r) (25)
where ~r and ~v are respectively the ECI position and velocity, and ω
~ E is the rotation
rate of the Earth. Assuming an atmosphere that rotates with the Earth, ω
~ E can be
expressed in terms of ECI components as:
 
0
 
ω
~E =  0  (26)
 
 
ωE
where ωE = 7.2921158553 × 10−5 rad/s [36].
Equation (24) can be inserted as part of the ~ad term in Eq. (22) to capture the
effect of aerodynamic (drag) perturbations. While Eq. (24) is a simple expression,
accurately modeling drag effects on satellites is a challenging prospect and an area
of ongoing research. The challenge in accurately modeling drag is largely due to the
difficulty in estimating the density (ζ). Atmospheric density has seasonal variations,
16
and is highly dependent on the current solar flux and Earth geomagnetic activity,
both of which are difficult to predict accurately. A discussion of many of the factors
at play when determining density can be found in ref. [25].
For the research presented in this thesis, the density is estimated using an
exponential model. This model only considers altitude variations, and assumes a
static (time-independent), axially symmetric atmosphere. In this case the density
can be estimated as:
γ − γref
ζ = ζref exp − (27)
Γ
where γ is the current altitude, γref is the reference altitude, ζref is the reference
density, and Γ the scale height. These parameters can be found in Table 1 [36].
Using Eq. (27) and Table 1, the density can be estimated as a function of γ, which
Table 1: Exponential Density Model

γ (km) γref (km) ζref (kg/m3 ) Γ (km)
200-250 200 2.789 × 10−10 38.5
−11
250-300 250 7.248 × 10 46.9
300-350 300 2.418 × 10−11 52.5
−12
350-400 350 9.158 × 10 56.4
400-450 400 3.725 × 10−12 59.4
−12
450-500 450 1.585 × 10 62.2
−13
500-600 500 6.967 × 10 65.8
600-700 600 1.454 × 10−13 79.0
−14
700-800 700 3.614 × 10 109.0
800-900 800 1.170 × 10−14 164.0
900 - 1,000 900 5.245 × 10−15 225.0
> 1,000 1000 3.019 × 10−15 268.0
can be approximated as γ = r − RE . Then, Eq. (24) can be used to estimate the

drag perturbation. The non-standard notation for density (ζ) and height (γ) in this
thesis is used to prevent any potential ambiguity with the bounded relative motion
parameter ρ or the orbit element h.
For satellite formations of identical spacecraft with the same attitude, the dif-
ferential drag (difference of drag force on the satellites) is essentially zero. In this
case, the spacecraft will have the same CD , S, and m, and the differences in ζ and
17
~vrel for the various satellites in the formation will be minuscule. This means that for
formations with identical design and attitude, the only effect drag will have on the
formation is the degradation of the reference (chief satellite) orbit.
However, there may be cases where there are formations with spacecraft whose
design is not identical, or which are unable to have the same attitude. In this case, it is
desirable to model the net effect of drag on the formation. The method implemented
in this thesis is not high-fidelity, but it does provide a reasoned estimate of differential
drag effects.
2.6 Mean Orbit Elements
In the two-body problem, the only classical orbit element that varies with time
is the anomaly (l or ν). The other five elements (a, e, i, h, and g) are constant.
However, in reality orbital perturbations cause these elements to vary with time.
These time-varying elements are referred to as the osculating elements.
Alternatively, mean orbit elements are orbit elements with the short period (on
the order of one orbit period) and long period (approximately one order higher than
one orbit period) effects averaged out. This property will be briefly demonstrated in
the following development [4]. The Delaunay variables are defined as:
D̄ = (l, g, h, L, GH) (28)
√ p
where L = µa, G = µa(1 − e2 ), and H = G cos i. The variables l, g, and h in this
case become generalized coordinates, and L, G, and H are their conjugate momenta.
The Delaunay elements are a set of canonical variables, meaning that they satisfy
Hamilton’s equations:
∂H ∂H ∂H
ġ = , l˙ = , ḣ = (29)
∂G ∂L ∂H
∂H ∂H ∂H
Ġ = − L̇ = − Ḣ = − (30)
∂g ∂l ∂h
18
The two-body Hamiltonian in terms of Delaunay variables is [4] [33]:
µ2
H0 = − (31)
2L2
This further demonstrates the fact that the mean anomaly l is the only non-constant
orbit element for the two-body problem, as all of the coordinates and the conjugate
momenta G and H do not appear in the two-body Hamiltonian. However, these results
are no longer valid in the presence of perturbations, as the Hamiltonian becomes:
H = H0 + H1 (32)
where H1 is the perturbing Hamiltonian and is the perturbation. In this case,

= −J2 , and the perturbed Hamiltonian is [4]:
µ4 Re2 a 3 h H 2 H2 i
H1 = 3 − 1 + 3 1 − cos θ (33)
4L6 r G2 G2
where θ is the argument of latitude (θ = ν + g). Note that H is now a function of all
of the coordinates and momenta except for h. Therefore, none of the coordinates or
momenta except for H will be constant.
Dirk Brouwer used the von Zeipel method in 1959 to average the Hamiltonian
in order to remove the periodic terms [37]. A simplified result of his original work
is presented here (Brouwer’s work included higher-order terms in the geopotential).
First the Hamiltonian is averaged with respect to l to remove the short-period terms,
then with respect to g to remove the long-period terms. If the distances are normalized
by the Earth radius, and time normalized by the mean motion of a satellite at one
Earth radius, the averaged Hamiltonians in terms of normalized units become [4]:
1
H̄0 = − (34)
2L̄2
1 L̄ H̄ 2
H̄1 = − 1 − 3 (35)
4L̄6 Ḡ Ḡ2
19
where the overhead bar on the Delaunay variables indicates that they are mean ele-
ments. The Hamiltonians are now a function of the momenta only (the coordinates
are ignorable). This means that the mean coordinate rates and the mean momenta
are constant, because the time derivatives of the momenta are zero (in accordance
with Hamilton’s equations), and the coordinate rates are a function of the momenta
only. This demonstrates the benefit to designing satellite formations in the mean
element space. In the osculating space, the coordinate rates and momenta are not
constant (observe Eq. 33), making it difficult to choose the correct initial values for a
formation. The caveat is that the use of mean elements requires a mean-to-osculating
transformation in order to get the instantaneous orbit elements. The actual trans-
formation from mean to osculating elements is obtained via a canonical coordinate
transformation using generating functions, a concept that will not be discussed in de-
tail here. The interested reader may find discussions of this topic in refs. [38] [39] [33].
Note that the development presented here is only to first-order, meaning that
Eq. (32) was truncated after the first expansion term. The perturbed Hamiltonian
could be expanded to any number of terms:
1
H = H0 + H1 + 2 H2 . . . (36)
2
However, the amount of algebra involving in expanding to higher terms quickly be-
comes nontrivial. In many cases a first-order expansion in J2 is adequate, as a
second-order J2 expansion is not necessary unless first-order expansions for higher-
order spherical harmonics are also needed [33]. The work done by Brouwer has been
extended by Gim and Alfriend, who derived the first-order mean-to-osculating trans-
formation for nonsingular elements [14] and equinoctial elements [40] (the D and D∗
matrices in the articles).
20
2.7 Small Satellites
Small Satellites are traditionally defined as the class of satellites with masses
below 500 kg. The SmallSat category can be broken in to four main sub-categories –
minisatellites, microsatellites, nanosatellites, and picosatellites [41]. These categories
are shown in Fig. 5. CubeSats are a specific type of SmallSat that are built using
standardized 10 cm × 10 cm × 10 cm cubes (1U). For example, a 6U CubeSat is a
satellite where the base structure is six 10-cm cubes. This standardization allows for
CubeSats to be developed and launched on a much faster timescale than traditional
satellites. CubeSats are typically nanosatellites, but can be in the microsatellite range.
For the purposes of this thesis, the term SmallSat can apply to any of the categories
of satellite in Fig. 5. The presented analysis will be focused mostly on satellites in
the microsatellite and nanosatellite range, with an emphasis on CubeSats.
Figure 5: Categories of Small Satellites
Satellites of the sizes illustrated in Fig. 5 tend to lend themselves towards

secondary payload launch options, as SmallSats are typically developed on a low-cost
budget. Additionally, there are now attractive options such as the Sherpa [21] for
SmallSats that provide additional ∆V capability after separation from the primary
launch vehicle, giving mission designers more control over the satellite’s final orbit.
In the realm of CubeSats, there are standardized options such as the canisterized
satellite dispenser (CSD) [22] that allow for relatively simple integration with launch
vehicles. These options available to SmallSats allow for quick development times and
less testing in comparison to larger satellites.
21
One drawback is that SmallSats typically have less control over their final or-
bit than a primary payload. This creates unique problems for the case of SmallSat
formation flying, since the spacecraft in the formation often have unfavorable initial
conditions, such as a large drift rate, that must be compensated for early in the mis-
sion. This is the primary motivation for the formation initialization algorithms that
will be developed in the Methodology chapter. Additional challenges for SmallSat
missions are low budgets and SWaP (size, weight, and power), which generally trans-
late to less capable subsystems. In effect, the same level of performance cannot be
expected of a SmallSat and a large satellite from a high-budget program. This leads
to challenges in the context of formation flying, as SmallSats will generally have lower
performance in terms of their navigation and propulsion subsystems as compared to
larger satellites.
2.8 State Transition Matrix
The state transition matrix (STM) is a matrix that maps the initial state X̄0
to the final state X̄f . Another way of describing it is the sensitivity matrix of the
current state to the initial conditions. The state transition matrix is defined as [33]:
∂ X̄(t)
Φ(t, t0 ) = (37)
∂ X̄(t0 )
So if the solution is known, the state transition matrix may be written in closed form.
The notation Φ(t, t0 ) defines the map from the initial state to the final state. For the
rest of this thesis, the initial time t0 is defined as zero, and the notation Φt is used.
22
For the Hill/Clohessy-Wiltshire model, the state transition matrix is [4]:
 
sin (nt)
4 − 3 cos (nt) 0 0 n
2
n
− 2 cosn(nt) 0
 
 −6nt + 6 sin (nt) 1 0 − n + 2 cosn(nt)
2 4 sin (nt)
− 3t 0
 
n

 
 sin (nt) 
 0 0 cos (nt) 0 0 n

Φt = 



 3n sin (nt) 0 0 cos (nt) 2 sin (nt) 0 
 
−6n + 6n cos (nt) −2 sin (nt) −3 + 4 cos (nt)
 
0 0 0 
 
0 0 −n sin (nt) 0 0 cos (nt)
(38)
Notice how the state transition matrix reduces to the identity matrix when t = 0.
State transition matrices only truly map the initial state to the final state for linear
systems. For nonlinear systems, it is maps the initial differential state to the final dif-
ferential state. This has applications in perturbation theory, as it can show how initial
trajectory errors will evolve over time [35] (departure from a reference trajectory).
Many of the available relative motion models are linear systems. In this case,
the state transition matrix can be used to determine the spacecraft state at some
future time based on the initial conditions:
X̄t = Φt X̄0 (39)
This makes the state transition matrix a powerful tool for trajectory propagation
and maneuver planning, as it is an analytical method that requires no numerical
integration. However, it must be noted that any linear relative motion model will be
an approximation of the true nonlinear dynamics. Therefore, it is important to be
sure that state transition matrices are only used where appropriate.
2.9 Maneuvering Techniques
There are a variety of propulsion methods available for spacecraft mission plan-
ners. Broadly, these can be separated into two categories – impulsive maneuvering
and continuous-thrust maneuvering. The spacecraft thrust can be modeled as impul-
23
sive if the thrust is high enough to impart the required velocity change in a nearly
instantaneous manner. In reality, all propulsion burns occur over some finite time
period. However, it is desirable if possible to model the burns as impulsive, as this
greatly reduces the complexity of the analysis.
In cases where the impulsive assumption is not valid, the spacecraft thrust may
be modeled as continuous. In this case, the thrust is modeled as an element of the
spacecraft dynamics, rather than a discontinuity in the spacecraft velocity [42] [43].
This point will be further emphasized in the following subsections.
2.9.1 Impulsive Maneuvering. In the impulsive maneuvering case, the ve-

locity change of the spacecraft (∆V ) is modeled as a discontinuity in the spacecraft
state. This is possible due to the assumption that the velocity change occurs in-
stantaneously (i.e. no time passes during the burn). Consider the state vector of a
spacecraft an instant before a maneuver:
 
~r0
X̄0 =   (40)
~v0
where ~r0 and ~v0 are respectively the initial position and velocity vectors of the space-
craft. The spacecraft state an instant after the maneuver can simply be modeled
as:  
~r0
X̄f =   (41)
~v0 + ∆~v
where ∆~v is the vector form of the velocity change (∆V = ||∆~v ||). The instantaneous
nature of the maneuver means that the initial (just before) and final (just after)
positions are the same, with jump discontinuities in the velocity.
Consider a general 2-impulse PCO reconfiguration maneuver. The horizontal

(y-z plane) component of this maneuver is shown in Fig. (6), with the first maneuver
at location “a” and the second maneuver at location “b”. Given the state transition
24
matrix Φ and the maneuver time tm , the required ∆V can be solved for as follows:
     
~rf ~
0 ~r0
X̄f =   =   + Φtm   (42)
~vf ∆~vb ~v0 + ∆~va
In this case, ~r and ~v are respectively the position and velocity of the spacecraft in
the Hill frame. The initial and final spacecraft states are known, and the ∆V s are
unknown. Therefore, Eq. (42) is a system of six equations with six unknowns that
can be solved for to find the required ∆~v components. For the impulsive maneuvers
in this thesis, the approach of Eq. (42) was modified by setting tm as a free variable
and optimizing for ∆V . This will be elaborated further in the Methodology chapter.
Figure 6: Impulsive Maneuver
2.9.2 Continuous-Thrust Maneuvering. In the case of continuous-thrust

maneuvering, the spacecraft thrust is not high enough to impart the required veloc-
ity change in an instantaneous manner. Examples of this are low-thrust propulsion
technologies such as ion drives or Hall thrusters that often require long burn times to
achieve the desired velocity change. Low-thrust technologies of this nature are often
25
desirable because their high specific impulse (Isp ) translates to low fuel consumption
as compared to chemical propulsion technologies.
In the continuous-thrust case, an analysis akin to Eq. (42) is no longer valid,

as the maneuver can no longer be treated as separate coasting phases with jump
discontinuities in the velocity between them. Instead, the thrust is added as a term
to the dynamics model:
X̄˙ = AX̄ + B~u (43)
In Eq. (43), A is the state-space matrix, B is the control matrix, and ~u is the control
acceleration. The matrix B is simply used to map the components of ~u onto the
correct equations of the state-space model.
As can be seen from Eq. (43), determining the required ~u throughout the
maneuver is not as straightforward as determining the required ∆~v in Eq. (42).
There are a variety of techniques that may be implemented to solve the continuous-
thrust maneuver problem. In low-thrust cases where the maneuver occurs over a
period of multiple orbits, the solution can be obtained by treating the control as a
perturbation [44].
Many of the references cited in the Introduction apply a variety of methods to

determine continuous-thrust solutions in the context of relative motion. The method
that was implemented for the research in this thesis is an LQR (Linear Quadratic
Regulation) technique derived by Cho and Park [16], and will be detailed further in
the Methodology chapter.
2.10 Propulsion Technologies
There are a variety of propulsion technologies that are available for spacecraft
missions. This section will overview two types of propulsion technologies – chemical
propulsion and electrical propulsion. Methods such as solar sails and nuclear propul-
sion will not be discussed. The comparison will be largely focused on performance
and design implications. The mechanics of how the different propulsion technologies
26
operate will not be discussed in detail. For a more complete discussion of electric
propulsion types, the reader is referred to [45]. Ref. [41] discusses many of the avail-
able propulsion technologies for small satellites. A large amount of the information
presented in this section was gleaned from these sources.
2.10.1 Chemical Propulsion. Chemical propulsion is the most commonly

implemented method of rocket propulsion, and the only method that is currently
feasible for launch vehicles, as alternative methods either do not generate enough
thrust or are not safe to operate in the atmosphere. Chemical propulsion methods
are also used on probe missions to the outer solar system, as they typically do not
require large amounts of power to operate.
Chemical propulsion methods have the advantage of high thrust as compared to

other propulsion technologies. It is for this reason that chemical propulsion methods
are often favored for large maneuvers, as the amount of time it would take for a
low-thrust method to impart the required velocity change is often prohibitive. The
disadvantage is that chemical propulsion methods have low specific impulse ratings
as compared to electrical propulsion systems. Specific impulses for SmallSat chemical
propulsion technologies range from 65 s to 250 s [41].
2.10.2 Electrical Propulsion. Electrical propulsion technologies have the

advantage of higher specific impulses than chemical propulsion methods. Typical
Is values for SmallSat technologies range from 700 s to 3000 s [41]. However, elec-
tric propulsion methods have much lower thrust ratings than chemical propulsion
methods, meaning that maneuver times are long. It is for this reason that electric
propulsion devices have typically been used for attitude control or orbit maintenance
on large satellites, rather than for performing large maneuvers.
For many small spacecraft, a high specific impulse is necessary to meet ∆V

requirements, as space for fuel storage is often limited. This may necessitate the
use of electric propulsion as the primary propulsion system. There are a variety of
27
electric propulsion technologies available, but this section will highlight three types –
ion propulsion, Hall thrusters, and electrospray thrusters.
2.10.2.1 Ion Propulsion. Ion propulsion uses plasma generation tech-

niques to ionize large fractions of the propellant, which is typically xenon gas. Ion
propulsion devices have the highest efficiency and specific impulse ratings of the pre-
sented electric propulsion technologies. Ion propulsion methods also require less power
as compared to Hall thrusters.
2.10.2.2 Hall Thrusters. Hall thrusters exploit a phenomenon known

as the Hall effect to generate plasma. Like ion propulsion devices, xenon gas is typ-
ically used as the propellant. Hall thrusters are less efficient and have lower specific
impulse ratings than ion propulsion devices. However, Hall thrusters have the highest
thrust of the presented electric propulsion methods, and are mechanically much sim-
pler than ion propulsion devices. The drawback for Hall thrusters is that they require
comparatively large amounts of power to operate. Hall thrusters for SmallSats have
power requirements ranging from 175-200 W, while power requirements for SmallSat
ion propulsion devices do not exceed 60 W [41].
2.10.2.3 Electrospray Thrusters. Electrospray propulsion devices elec-

trostatically extract and accelerate ions from a conductive salt with a negligible vapor
pressure. These devices have somewhat lower thrust ratings than ion propulsion meth-
ods, and comparable specific impulse ratings to Hall thrusters. While this may make
them seem unappealing from a performance standpoint, electrospray devices have
many advantages as compared to the other electric propulsion technologies. One of
the biggest advantages for electrospray devices is that the propellant does not need to
be pressurized for storage. Additionally, power requirements to operate electrospray
devices are comparatively low ( < 15 W) [41]. These design advantages make elec-
trospray devices an attractive option for SmallSat applications, where storage space
and power are often limited.
28
III. Methodology
This chapter will outline the methods used to generate the results presented in this
thesis. The GCO initial conditions are derived in terms of equinoctial elements,
followed by a description of the methods for computing fuel-optimal trajectories.
Transformations between ECI and Hill frame components are detailed. Methods are
developed for negating relative drift in the context of formation initialization. The full
nonlinear equations of motion with perturbations are provided. Finally, a top-level
overview of the algorithms is presented to show how the various pieces are connected.
3.1 Formation Design
In order to establish or reconfigure a spacecraft formation, it is necessary to first

determine the initial conditions that will result in the desired formation geometry. It
is advantageous to express these initial conditions in terms of orbit element differences,
as the Hill frame components are all fast variables, making it difficult to determine
the correct value at a given point in time. It is worth reiterating that the initial
conditions should be used as the mean element differences when designing formations
in the presence of perturbations. The PCO and GCO initial conditions in terms
of nonsingular elements were derived by Vaddi, Alfriend, Vadali, and Sengupta [11]
Their work was extended by Johnson to derive the PCO initial conditions in terms of
equinoctial element differences [28].
The following development will build on the previous work by deriving the
GCO initial conditions in terms of equinoctial elements so that GCOs may be used in
equatorial orbits. The parameters for bounded relative motion in terms of nonsingular
orbit element differences are provided in Eqs. (44-48) [28].
ρx = ac δe (44)
ρy = ac (δλ0 + δh cos (ic )) (45)

p
ρz = ac δi2 + δh2 sin2 ic (46)
29
π
αx = λc0 − atan2(δq2 , δq1 ) − (47)
2
αz = λc0 + atan2(−δh sin ic , δi) (48)
where the subscript c refers to the chief satellite and the subscript 0 means the value
at the initial time. The parameters defined in Eqs. (44 - 48) are the same parameters
that were previously given in terms of Hill frame components (11 - 15). The GCO
constraints are defined as [11] [4]:
δa = 0 (49)
ρy = 0 (50)
ρz
ρx = √ (51)
3
αx = α + λc0 (52)
αz = α + λc0 (53)
where α is the deputy phase angle at the time of the chief satellite’s equator crossing.
These constraints may then be written in terms of the nonsingular element differences:
δa = 0 (54)
π
−atan2(δq2 , δq1 ) − =α (55)
2
atan2(−δh sin ic , δi) = α (56)
δλ0 + δh cos ic = 0 (57)

ρz
δe = √ (58)
ac 3
p ρz
δi2 + δh2 sin2 ic = (59)
ac
Assuming a circular chief satellite orbit, the eccentricity of the deputy is equal to the
eccentricity difference:
ed = δe (60)
30
Therefore:
q
δe = δq12 + δq22 (61)
The assumption of a circular chief orbit is valid, because GCO (and PCO) forma-
tions are defined in terms of the HCW solution. Therefore, they only truly exist
for circular chief orbits. However, the GCO (and PCO) initial conditions will re-
sult in stable motion for elliptical chief orbits (but the resulting formation will not
have the unique properties of a PCO or GCO). Using Eq. (61), and the identity
π
−atan2(δq2 , δq1 ) − 2
= atan2(−δq1 , −δq2 ), Eq. (55) and Eq. (58) may be re-written:
atan2(−δq1 , −δq2 ) = α (62)
ρz
q
δq12 + δq22 = √ (63)
ac 3
Equations (62) and (63), along with Eqs. (56), (57), and (59), are the equations of
constraint for a GCO in terms of nonsingular element differences. Equation (54) is
the previously discussed no-drift condition, and should be replaced by an appropriate
expression in the presence of perturbations [46] [4] [47].
In order to avoid the equatorial orbit singularity associated with nonsingular

elements, the equations of constraint may be written in terms of differences in the
equinoctial elements. The transformations from nonsingular elements to equinoctial
elements are:
λ=Λ−h (64)
δλ = δΛ − δh (65)
q
i = 2atan( p21 + p22 ) (66)
2(p1 δp1 + p2 δp2 )

δi = p (67)
p21 + p22 (1 + p21 + p22 )
h = atan2(p2 , p1 ) (68)
p1 δp2 − p2 δp1
δh = (69)
p21 + p22
31
1
δq1 = p 2 δh(q̃ 2 p1 − q̃1 p2 ) + δ q̃1 p1 + δ q̃2 p2 (70)
p1 + p22
−1
δq2 = p 2 δh(q̃ 1 p1 + q̃2 p2 ) + δ q̃1 p2 − δ q̃2 p1 (71)
p1 + p22
Assuming zero chief eccentricity:
1
δq1 = p 2 (δ q̃1 p1 + δ q̃2 p2 ) (72)
p1 + p22
−1
δq2 = p 2 (δ q̃1 p2 − δ q̃2 p1 ) (73)
p1 + p22
Using these transformations on Eqs. (56-57, 59, 62-63), and introducing a new phase
angle αI = α − hc , yields the equations of constraint in terms of equinoctial elements
and equinoctial element differences. αI is the deputy phase angle when the chief
satellite crosses the 1-direction of the ECI (Earth-Centered Inertial) frame, known as
the first point of Aries.
p1c δp2 − p2c δp1

q
δΛ0 = 2 2
[1 − cos 2atan( p21c + p22c ) ] (74)
p1c + p2c
atan2(−(p1c δp2 − p2c δp1 ), p1c δp1 + p2c δp2 ) = αI + atan2(p2c , p1c ) (75)
s
1 4(p1c δp1 + p2c δp1 )2 (p1c δp2 − p2c δp1 )2 ρz
q
2 2 2
p 2 2 2
+ 2 2
sin (2atan( p 1c + p 2c )) =
p21c + p22c (1 + p1c + p2c ) p1c + p2c ac
(76)

atan2 − (δ q̃1 p1c + δ q̃2 p2c ), −(−δ q̃1 p2c + δ q̃2 p1c ) = αI + atan2(p2c , p1c ) (77)
Since the chief satellite’s orbit is assumed to be circular, δe = ed , so the final equation
of constraint is:
ρz
q
δ q̃12 + δ q̃22 = √ (78)
ac 3
After some algebra, these equations can be solved to yield the GCO initial conditions
in equinoctial element differences:
δa = 0 (79)
32
ρz
δ q̃1 = − √ sin αI (80)
ac 3
ρz
δ q̃2 = − √ cos αI (81)
ac 3
ρz
δp1 = (1 + p21c + p22c ) cos αI (82)
2ac
ρz
δp2 = − (1 + p21c + p22c ) sin αI (83)
2ac
2(p1c δp2 − p2c δp1 )
δΛ0 = (84)
1 + p21c + p22c
However, it is desirable to express these equations in terms of a more intuitive physical
property. Using the knowledge that the radius of a GCO is equal to 2ρx , a variable
ρ = 2ρx = √2 ρz may be defined. Then, the GCO initial conditions terms of ρ are:
3
ρ
δ q̃1 = − sin αI (85)
2ac
ρ
δ q̃2 = − cos αI (86)
2ac
√
3ρ
δp1 = (1 + p21c + p22c ) cos αI (87)
4 ac
√
3ρ
δp2 = − (1 + p21c + p22c ) sin αI (88)
4 ac
Equations (85-88), along with Eqs. (79) and (84), will yield the necessary initial
conditions in equinoctial element differences to yield a GCO of desired radius ρ and
angle αI . As mentioned previously, Eq. (79) should be replaced by an appropriate
term in the presence of perturbations. The δa to negate drift in presence of J2
perturbations given by Eq. (89) [4]. This constraint is also valid in the two-body
problem, as it reduces to δa = 0 when J2 vanishes.
R 2 3η + 4
e c
δa = 0.5J2 ac [(1 − 3 cos2 ic )δη − (ηc sin 2ic )δi] (89)
ac ηc5
33
where
p
ηc = 1 − e2c (90)
ec δe
δη = − (91)
ηc
The equations for the GCO initial conditions in terms of equinoctial element
differences are summarized in Table 2, along with the initial conditions for a GCO in
terms of nonsingular elements [11], and the PCO initial conditions in nonsingular and
equinoctial elements [11] [28]. For all of the cases in Table 2, an additional equation of
Table 2: Summary of PCO and GCO Initial Conditions

PCO GCO
√
δλ = aρc sin α
tan i
δλ = 23 aρc sin α
tan i
√
δi = aρc cos α δi = 23 aρc cos α
ns
¯ δq1 = − 2aρ c sin α δq1 = − 2aρ c sin α
δq2 = − 2aρ c cos α δq2 = − 2aρ c cos α
√
δh = − aρc sin α
sin i δh = − 23 aρc sin α
sin i
δ q̃1 = − 2aρ c sin αI δ q̃1 = − 2aρ c sin αI
δ q̃2 = − 2aρ c cos αI δ q̃2 = − 2aρ c cos αI
√
ρ 2 2
¯ δp1 = 2ac (1 + p1c + p2c ) cos αI
eq δp1 = 43 aρc (1 + p21c + p22c ) cos αI
√
δp1 = 2aρ c (1 + p21c + p22c ) sin αI δp2 = − 43 aρc (1 + p21c + p22c ) sin αI
δp2 −p2c δp1 )
δΛ0 = 2(p1c1+p 2 +p2 δΛ0 = 2(p1c1+pδp2 −p2c δp1 )
2 +p2
1c 2c 1c 2c
constraint (Eq. (89)) is needed to negate relative drift. As was discussed previously,
the design parameter ρ defines the radius of a GCO. For PCOs, the parameter ρ
defines the radius of the projected circle in the horizontal plane. The parameters α
and αI have the same meaning as discussed previously.
3.2 Impulsive Formation Control
This section will discuss the methods used for the impulsive thrust algorithm.
First, the method for formation initialization and reconfiguration will be discussed,
followed by a discussion of the formation maintenance algorithm.
34
3.2.1 Impulsive Maneuvering Algorithm. The method used to solve for
impulsive thrust maneuvers is a modification of Eq. (42). Rather than specifying the
transfer time and solving for the ∆V s, the transfer time was set as a free optimization
variable, subject to the inequality constraint:
∆tm < ∆tmax (92)
where ∆tmax is the maximum allowable time for the maneuver. A direct optimization
routine using MATLAB’s fmincon function was used to search transfer time values to
find the fuel-optimal (∆V -optimal) maneuver. It is desirable to analyze cases involv-
ing more than two burns, so Eq. (42) was programmed recursively for N maneuvers
with N − 1 coasting arcs:
         
~rf ~
0 ~
0 ~
0 ~r
X̄f =   =   + Φ∆tm   + Φ∆tm   + Φ∆tw  0  . . . (93)
2 1
~vf ∆~vc ∆~vb ∆~va ~v0
where Φ∆tmi is the STM for the ith coasting arc, and ∆~va , ∆~vb , ∆~vc and so on are the
N burns. Φ∆tw is the STM for an initial “wait” coast, since it may not be optimal to
begin maneuvering immediately. The equations of constraint are:
X̄f = X̄fd (94)
(N −1)
X
∆tm = ∆tw + ∆tmi < ∆tmax (95)
i=1
35
where X̄fd is the desired final state. ∆tw is the time duration of the wait period, and
∆tmi is the duration of the ith coast arc. The optimization state vector is:
 
∆~v1
 .. 
 
 . 
 
 ∆~vN 
 
 
x̄opt =  ∆tm1  (96)
 
 .. 
 
 . 
 
 
∆tmN −1 
 
∆tw
with the cost function defined as:
J = x̄Topt M x̄opt (97)
..
 
I3N .
M =  (98)
. . . 0N
where I3N is a 3N × 3N identity matrix. Since there are 3 components for each ∆~v ,
M is a matrix that maps only the ∆~v components to the quadratic cost function. The
cost function is minimized, subject to Eqs. (94) and (95), resulting in a minimum-fuel
impulsive maneuver.
This algorithm has the advantage of being general. Any number of N impulses
could be defined, allowing a wide variety of formation reconfiguration and initialization
maneuvers to be analyzed. Additionally, any state transition matrix can be used with
this method (as long as the state vector is defined with respect to the model). The
disadvantage is that it is not as computationally efficient as an analytic method. The
optimization state vector has 4N components, so increasing the number of maneuvers
rapidly increases the number of optimization variables. This means that simulating
cases with large N values is computationally prohibitive. The state transition matrix
used for the research in this thesis is the STM for the Schweighart-Sedwick model [26],
36
in order to capture the effects of J2 in the targeting dynamics. It is provided in
Appendix A.
3.2.2 Impulsive Formation Maintenance Algorithm. The method used for

impulsive formation initialization and reconfiguration maneuvers has the advantage
of allowing a wide variety of scenarios to be analyzed. While the method could con-
ceivably be used for formation maintenance, it was decided to implement a more
computationally efficient approach, since formation maintenance analyses usually re-
quire long run times.
The approach used for formation maintenance in the impulsive-thrust case is a

discrete-time linear quadratic regulator (DLQR). The method implemented for this
thesis is a modified version of the DLQR approach outlined in ref. [4]. A DLQR is a
discretized version of the continuous-time linear quadratic regulator (LQR). DLQRs
(and LQRs) are an analytic control approach that uses linearized dynamics (such as
the HCW or Schweighart-Sedwick equations) and a quadratic cost function. The cost
function for the DLQR is given by Eq. (99).
n
X
J= X̄(i)T QX̄(i) + ~u(i)T R~u(i) (99)
i=0
where n is the number of time steps, and X̄(i) and ~u(i) are the state and control
vectors at the ith time instant. Q and R are respectively the state and control weight
matrices. The state and control vectors are subject to the discrete-time state equation:
X̄(i + 1) = Ad X̄(i) + Bd~u(i) (100)
Note that even though Eq. (100) takes a similar form to Eq. (43), the quantities
are fundamentally different. The left side of Eq. (100) is the state at a future time
instant, not the state derivative. This means that Ad is not the state-space matrix,
37
but rather the state transition matrix (assuming impulsive control):
Ad = Φ∆t (101)
where ∆t is the time between instants (the step size). Bd then takes the form:
Bd = Ad B = Φ∆t B (102)
where B is the same control matrix that was discussed previously. The approximate
control law is given by ref. [4]:

~u(i) = −K X̄(i) − X̄r (i) (103)
where K = (Bd SBd + R)−1 BdT SAd is the steady-state Kalman gain matrix, and X̄r is
the reference (desired) trajectory. The control ~u when using a DLQR is actually the
∆~v , rather than a control acceleration vector as in the continuous-time case. S is the
solution to the discrete-time algebraic Riccati equation:
ATd SAd − S − ATd SBd (BdT SBd + R)−1 BdT SAd + Q = 0 (104)
Numerical solvers of the discrete-time algebraic Riccati equation (Eq. (104)) are
1
standard in many software packages, such as MATLAB’s dare function. The
Schweighart-Sedwick STM was used as the Ad matrix, and X̄r was determined by
sampling the reference trajectory of the desired formation at the instant of the im-
pulse application.
The DLQR is an optimal control technique that has the advantage of being
a computationally efficient analytic method for formation maintenance. Equation
(104) must be solved numerically, but the resulting gain matrix is then constant as
1
MATLAB R is a registered trademark of MathWorks, Inc.
38
long as the time intervals between impulses remains the same. This means that the
discrete-time algebraic Riccati equation only needs to be solved once.
3.3 Continuous-Thrust Formation Control
This section discusses the methodology that was implemented for computing
low-thrust solutions in relative motion. The method was taken from Cho and Park
[16], with some modifications that were made to constrain the solutions. The reader
is encouraged to read ref. [16] for a more comprehensive discussion of the method.
The continuous-thrust control is an LQR (linear-quadratic regulator) approach,

assuming three-axis control and a specified transfer time. The method can be imple-
mented for a variety of dynamic models, as long as the model uses linearized dynamics
such that the elements of the state-space matrix are related by A1 − AT1 = Ȧ2 .
 
03 I3
A=  (105)
A1 A2
For the research presented in this thesis, the Schweighart-Sedwick model [26] was
used to capture the effects of J2 in the thrust profile calculation. Given the state-
space matrix A and state transition matrix Φ for the desired dynamic model, the
fuel-optimal control can be calculated using the following equations:
Z t
S(t) = ΦTA ΦA dτ (106)
0
K = Φ−1 −1
f X̄f − Φ0 X̄0 (107)
C = ΦTA Φ̇A − (ΦTA Φ̇A )T − ΦTA A2 ΦA (108)
where
X̄t = Φt X̄0 (109)
 
ΦA
Φ=  (110)
Φ̇A
39
 
x
 
y 
 
 
 
z 
X̄ = 
 
 (111)
ẋ
 
 
ẏ 
 
ż
X̄˙ = AX̄ + B~u (112)

 
03
B=  (113)
I3
the cost function is: Z t

1
J= uT u dτ (114)
2 0
and the optimal control is:

~u(t) = ΦA Sf−1 CK (115)
The matrices C and K are constant, and Sf is the S matrix evaluated at the
final time. Therefore, the optimal control ~u(t) throughout the transfer is a function
of ΦA (t) only.
The method of Cho and Park was modified by placing constraints on the achiev-
able trajectories. This can be done by bounding the control acceleration by the max-
imum thrust and initial mass of the spacecraft (U = ||~u||).
Fmax
Umax = (116)
m0
The maximum thrust is a constant, and the mass of the spacecraft can only decrease
throughout the maneuver. Therefore, the constraint in Eq. (116) is adequate for
bounding the control throughout any maneuvers. For long maneuvers with signifi-
cant changes in spacecraft mass, it may be desirable to update Umax throughout the
maneuver.
40
3.4 Coordinate Transformations
In order to be able to effectively implement formation flying algorithms, it is

necessary to be able to accurately transform state vectors between inertial and relative
motion frames. This section will provide the coordinate transformations between ECI
and Hill frame components, and the conversion of equinoctial elements to an ECI
state vector.
3.4.1 ECI to Hill Frame Transformation. The relative motion between two
satellites in general elliptic orbits can be expressed in the Hill frame as [48]:
δ~rT ~rc
x= (117)
rc

δ~r T ~ c × ~rc
H
y= (118)
~ c × ~rc ||
||H
~c
δ~rT H
z= (119)
Hc
~ is the angular momentum vector:
where H
~ = ~r × ~v
H (120)
~ δ~r and δ~v are respectively the relative displacement and velocity:
and H = ||H||.
δ~r = ~rd − ~rc (121)
δ~v = ~vd − ~vc (122)
Note that δ~r and δ~v are simply the difference in the ECI position and velocity be-
tween the deputy and chief, not Hill frame components. Equations (117-119) can be
differentiated with respect to time to yield the Hill frame velocity components:
δ~v T ~rc + δ~rT ~vc (δ~rT ~rc )(~rcT ~vc )

ẋ = − (123)
rc rc3
41
~ c × ~rc ) + δ~rT (H
δ~v T (H ~˙ c × ~rc + H ~ c × ~vc )
ẏ =
||H~ c × ~rc ||
δ~rT (H~ c × ~rc )(H

~ c × ~rc )T (H ~˙ c × ~rc + H~ c × ~vc )
− (124)
||H~ c × ~rc ||3
δ~v T H ~˙ c δ~rT H
~ c + δ~rT H ~ c (H
~ TH~˙ c )
c
ż = − (125)
Hc Hc3
~˙ c = ~rc × ~v˙ c = ~rc × ~ac . This means that determining the Hill frame velocities
where H
requires knowledge of the chief satellite’s acceleration vector, which can be obtained
from Eq. (187). Equations (117-119) and (123-124) are the exact transformation from
ECI to Hill frame components, and they are valid in the presence of perturbations.
3.4.2 Hill Frame to ECI Transformation. The rotation matrix from the Hill
to ECI frame can be expressed as:
h i
T = r̂c (Ĥc × r̂c ) Ĥc (126)
where r̂c and Ĥc are respectively the unit vectors in the direction of the chief ECI
position vector and angular momentum vector. This can be used to rotate the control
vector from the Hill to ECI frame:
~uECI = T ~u (127)
Equation (127) can be inserted into Eq. (186) to find the total ECI acceleration of a
deputy satellite during a low-thrust maneuver.
A similar process may be used to find the deputy ECI position and velocity
from its Hill frame components. Figure 7 shows the orientation of the deputy ECI
velocity vector in space.
42
Figure 7: Spacecraft Velocity Relative to ECI Frame
From the transport theorem, the velocity of the deputy with respect to the ECI
frame can be expressed as [49]:
~vd/ECI = Hill~r˙d/Hill + ω
~ Hill/ECI × ~rd/Hill + ECI ~r˙Hill/ECI (128)
where Hill
~r˙d/Hill is the time derivative in the Hill frame of the position of the deputy
with respect to the Hill frame, ω
~ Hill/ECI is the angular velocity of the Hill frame with
respect to the ECI frame, ~rd/Hill is the position of the deputy with respect to the Hill
frame, and ECI
~r˙Hill/ECI is the time derivative in the ECI frame of the position of the
Hill frame with respect to the ECI frame.
Equation (128) is a coordinate-free representation of the deputy velocity with

respect to the ECI frame. If ~vd/ECI and ECI
~r˙Hill/ECI are expressed in terms of ECI
components, then ~vd/ECI = ~vd and ECI ~r˙Hill/ECI = ~vc , where ~v is a vector of the velocity
components of the ECI state. Additionally, if ~rd/Hill and Hill
~r˙d/Hill are expressed in
terms of Hill frame components, then ~rd/Hill = ~rdHill and Hill
~r˙d/Hill = ~vdHill (where
~rdHill and ~vdHill are respectively the position and velocity components of the Hill-frame
state). If ω
~ Hill/ECI is also expressed in terms of Hill frame components, then Eq. (128)
can be modified:

~ Hill/ECI × ~rdHill + ~vc
~vd = T~vdHill + T ω (129)
43
The deputy ECI position is simply:
~rd = ~rc + T ~rdHill (130)
Equations (129) and (130) can be used to find the deputy ECI position and velocity,
based on the deputy Hill frame state and the chief ECI position and velocity. The
only remaining step is to determine ω
~ Hill/ECI . In the Hill frame, it can be expressed
as:  
ω
 r
ω
~ Hill/ECI =0 (131)
 
 
ωn
where ωr and ωn are respectively the radial and normal components of the angular
velocity. In the two-body problem, ωr is zero, and ωn is the mean motion. However,
J2 perturbations cause a radial component to exist for the angular velocity, due to
nodal regression. ωr and ωn can be expressed in terms of equinoctial elements as [40]:
√
2 µσ3 τ2
ωr = −3J2 Re2 7/2 (1 + q̃1c cos (Ψc ) + q̃2c sin (Ψc ))3 (132)
ac η 7 σ22
r
µ
ωn = (1 + q̃1c cos (Ψc ) + q̃2c sin (Ψc )) (133)
p
Ψ = ν + g + h is the equinoctial element analogous to the true anomaly. Conversions
between Λ and Ψ can be done with a modified version of Kepler’s equation, which
is discussed in the following section. The remainder of the parameters in Eqs. (132)
and (133) are [40]:
τ2 = p1c sin (Ψc ) − p2c cos (Ψ) (134)
q
2 2
η = 1 − q̃1c − q̃2c (135)
σ2 = 1 + p21c + p22c (136)
σ3 = 1 − p21c − p22c (137)
44
2 2

p = ac 1 − q̃1c − q̃2c (138)
The osculating equinoctial elements should be used to calculate ωr and ωn , and the
parameters τ2 , η, σ2 , σ3 , and p. With ωn and ωr defined, Eq. (129) may be rewritten
in matrix form:
       
Ẋ Ẋ ẋ −ωn y
 d  c    
~vd = [vd ] =  Ẏd  =  Ẏc  + T ẏ  + ωn x − ωr z  (139)
       
       
Żd Żc ż ωr y
The process in this section can also be used to rotate a ∆~v into the ECI frame
to apply it to the deputy ECI state. This is done by simply updating the deputy Hill
frame velocity to include the ∆~v before using Eq. (139).
3.4.3 Determining ECI Position and Velocity from Equinoctial Elements.

The ECI position and velocity of a satellite can be expressed in terms of its equinoctial
elements as [31]:
~r = X1 f~ + Y1~g (140)
~v = Ẋ1 f~ + Ẏ1~g (141)
where  
1− p22 + p21
1  
f~ = (142)
 
2p1 p2
1 + p21 + p22
 
 
−2p2
 
2p1 p2
1 
2 2

~g = 1 + p2 − p1  (143)
 
1 + p21 + p22  
2p1
and
X1 = a (1 − q̃22 β) cos (F ) + q̃2 q̃1 β sin (F ) − q̃1

(144)
Y1 = a (1 − q̃12 β) sin (F ) + q̃2 q̃1 β cos (F ) − q̃2

(145)
45
a2 n
Ẋ1 = (q̃2 q̃1 β cos (F ) − (1 − q̃22 β) sin (F )) (146)
r
a2 n
Ẏ1 = (−q̃2 q̃1 β sin (F ) + (1 − q̃12 β) cos (F ) (147)
r
1
β= p (148)
1 + 1 − q̃12 − q̃22
r
µ
n= (149)
a3
where F is the eccentric longitude, and is found by solving the modified Kepler’s
equation:
F − q̃1 sin (F ) + q̃2 cos (F ) = Λ (150)
The modified Kepler’s equation can be solved with the following algorithm:
while ∆F > tol
F − q̃1 sin (F ) + q̃2 cos (F ) − Λ

Fnew = F −
1 − q̃1 cos (F ) − q̃2 sin (F )
∆F = Fnew − F
F = Fnew
where tol is a specified numerical tolerance. The numerical tolerance is necessary

because Eq. (150) is a transcendental equation (as is the classical Kepler’s equation),
meaning that it cannot be solved explicitly. An initial guess of F = Λ may be used.
Once F has been found, it can also be used to find Ψ:
Ψ = atan2(ψ1 , ψ2 ) (151)
ψ1 = ((1 + η)(η sin (F ) − q̃2 ) + q̃2 (q̃1 cos (F ) + q̃2 sin (F )) (152)
ψ2 = (1 + η)(η cos (F ) − q̃1 ) + q̃1 (q̃1 cos (F ) + q̃2 sin (F )) (153)
46
3.5 Initial Conditions for Formation Initialization
In order to properly address the formation initialization problem for SmallSats,

the initial state of the spacecraft should be chosen to simulate a realistic dispersal
scenario. Figure 8, taken from the CSD data sheet [22], provides the ejection velocity
as a function of spacecraft mass.
Figure 8: CSD Ejection Velocity
In order to be able to establish the desired formation, the deputy spacecraft

must be able to cope with potentially unfavorable scenarios when it is deployed. A
dispersal with velocity in the ŷ direction would result in a large drift rate due to
the secular term (Eq. 5). Dispersals in the z-direction would result in cross-track
oscillatory motion but no drift. Similarly, a dispersal in the x-direction would result
in a 2:1 bounded relative motion ellipse going through the original dispersal point.
Both of these latter cases have the advantage of keeping the deputy relatively close
to the dispersal vehicle, but risk of collisions would be high if the formation was left
uncontrolled. This seems to indicate that in cases where the dispersal vehicle and the
chief satellite are the same, the ideal case would be aligning the dispersal mechanism
(such as the CSD) mostly along the x- or z-direction, with a small component in
the y-direction. This would have the advantage of keeping the drift rate small while
ensuring that no collision would occur.
47
If the chief and the deputy are both being dispersed from a carrier vehicle (as
opposed to a mothership chief satellite deploying its deputy satellites), the drift rate
between the two can be mitigated by deploying the satellites at one orbit intervals.
If both satellites are deployed along the y-axis at one-orbit intervals, the relative
velocity between the satellites will be close to zero. This result is demonstrated in the
following development. The HCW STM evaluated at t = P (corresponding to one
orbit period) is:  
1 0 0 0 0 0
 
−6nP 1 0 0 −3P 0
 
 
 
 0 0 1 0 0 0
Φ2π =


 (154)
 0 0 0 1 0 0
 
 
 0 0 0 0 1 0
 
0 0 0 0 0 1
If this is used to propagate a dispersal along the y-axis, the result is:
    
1 0 0 0 0 0
0 0
    
−6nP 1 0 0 −3P 0  0  −3P vdisp 
    
    
    
 0 0 1 0 0 0   0   0 
Φ2π X̄0 = 



=
 

 (155)
 0 0 0 1 0 0  0   0 
    
    
 0 0 0 0 1 0 vdisp   vdisp 
    
0 0 0 0 0 1 0 0
where vdisp is the dispersal velocity. If the satellites are of equal design they should
have approximately the same dispersal velocity. If a second satellite, defined as the
48
chief, is deployed one orbit later, its initial Hill frame state is given by:
 
0
 
 0 
 
 
 
 0 
X̄c0 = 
 
 (156)
 0 
 
 
vdisp 
 
0
It is desirable to express the relative motion of the deputy with respect to the
chief, rather than the dispersal vehicle. The Hill frame state of the deputy (in a new
Hill frame centered on the chief) can be found by determining the ECI state vector of
each spacecraft, and then using them to determine the state vector of the deputy in
the new Hill frame. The ECI position and velocity of the chief and deputy spacecraft
can be found using Eqs. (129) and (130).
~rc = ~rv + T ~rcHill = ~rv (157)
~vc = T~vcHill + T (~ω × ~rcHill ) + ~vv = T~vcHill + ~vv (158)
~rd = ~rv + T ~rdHill (159)
~vd = T~vdHill + T (~ω × ~rdHill ) + ~vv (160)
where the subscript v refers to the dispersal vehicle. ~r and ~v vectors with a “Hill”
subscript indicate the position and velocity components of the Hill frame state vector.
The position and velocity differences can be found using Eqs. (121) and (122):
δ~r = ~rd − ~rc = T ~rdHill (161)
δ~v = ~vdECI − ~vcECI = T (~vdHill − ~vcHill ) + T (~ω × ~rdHill ) (162)
49
However, ~vdHill = ~vcHill in this case, therefore:
δ~v = T (~ω × ~rdHill ) (163)
T is given by Eq. (126), and ~rdHill is known from the deputy Hill frame state vector.
If two-body motion is assumed, the angular velocity vector in the Hill frame is:
 
0
 
ω
~ = 0 (164)
 
 
n
where n is the mean motion. Now, δ~r and δ~v can be expressed as:
 
0
h i 
δ~r = r̂v Ĥv × r̂v Ĥv −3P vdisp  = −3P vdisp Ĥv × r̂v (165)
 
 
0
  
0 −n 0 0
h i  
δ~v = r̂v Ĥv × r̂v Ĥv n 0 0 −3P vdisp 
  
  
0 0 0 0
 
3P nvdisp
h i 
= r̂v Ĥv × r̂v Ĥv   = 3P nvdisp r̂v (166)
 
0
 
0
Now that the position and velocity differences are known, Eqs. (117-119) and
(123-125) can be used to find the Hill frame components. The transformation will
give the state of the deputy in a Hill frame centered on the chief satellite (where
previously it was centered on the dispersal vehicle). The development will start with
50
the x-component:
T
δ~rT ~rc −3P vdisp Ĥv × r̂v ~rc
x= = (167)
rc rc
However, ~rc = ~rv because the transformation is made at the instant the chief is ejected
from the dispersal vehicle. Eq. (167) can be rewritten as:
T
x = −3P vdisp Ĥv × r̂c r̂c = 0 (168)
Now, the y-component:
T
δ~rT H~ c × ~rc −3P vdisp Ĥv × r̂c ~ c × ~rc
H
y= = (169)
~ c × ~rc ||
||H ~ c × ~rc ||
||H
The cross product magnitude of two arbitrary vectors is ||~a ×~b|| = ab sin Θ, where Θ is
~ c and ~rc are orthogonal, ||H
the angle between the two vectors. Since H ~ c ×~rc || = Hc rc .
Therefore:
T
y = −3P vdisp Ĥv × r̂c Ĥc × r̂c (170)
The inner product of two vectors is ~aT~b = ab cos Θ. Assuming that the dispersal
vehicle and the chief satellite lie in the same orbit plane (which should be the case
T
for a y-axis dispersal), Ĥc = Ĥv . Therefore, Ĥv × r̂c Ĥc × r̂c = 1, because

the angle between the unit vectors Ĥv × r̂c and Ĥc × r̂c is zero. This gives the
expression for the y-component:
y = −3P vdisp (171)
The z component is:
~c
δ~rT H T
z= = −3P vdisp Ĥv × r̂c Ĥc = 0 (172)
Hc
51
The ẋ component is given by:
δ~v T ~rc + δ~rT ~vc (δ~rT ~rc )(~rcT ~vc )

ẋ = −
rc rc3
T
T
T
(3P nvdisp r̂c ) ~rc − 3P vdisp Ĥv × r̂c ~vc −3P vdisp Ĥv × r̂c ~rc (~rcT ~vc )
= −
rc rc3
T
T
(3P nvdisp r̂c )T ~rc − 3P vdisp Ĥv × r̂c ~vc 3P vdisp Ĥv × r̂c r̂c (~rcT ~vc )
= +
rc rc2
(173)
T
The 2nd term goes to zero, because Ĥv × r̂c r̂c = 0.
T
(3P nvdisp r̂c )T ~rc − 3P vdisp Ĥv × r̂c ~vc
ẋ =
rc
T
3P vdisp Ĥv × r̂c ~vc
= 3P nvdisp r̂cT r̂c −
rc
T
3P vdisp ~vc × Ĥv r̂c
= 3P nvdisp − (174)
rc
rc ×~vc
~
Ĥv = Ĥc = ||~
rc ×~vc ||
, therefore:
T
rc ×~vc
~
3P vdisp ~vc × ||~
rc ×~vc ||
r̂c
ẋ = 3P nvdisp −
rc
3P vdisp (~vc × ~rc × ~vc )T r̂c
= 3P nvdisp −
rc ||~rc × ~vc ||
T
3P vdisp ~rc ~vcT ~vc − ~vc ~vcT ~rc r̂c
= 3P nvdisp − (175)
rc ||~rc × ~vc ||
52
If the chief orbit is circular, ~vc and ~rc are orthogonal, so ~vcT ~rc = 0:
T
3P vdisp ~rc ~vcT ~vc r̂c
ẋ = 3P nvdisp −
rc ||~rc × ~vc ||
3P vdisp vc2~rcT r̂c
= 3P nvdisp −
rc ||~rc × ~vc ||
3P vdisp vc2 r̂cT r̂c
= 3P nvdisp −
||~rc × ~vc ||
3P vdisp vc2
= 3P nvdisp − (176)
||~rc × ~vc ||
Since ~rc and ~vc are orthogonal, ||~rc × ~vc || = rc vc :
3P vdisp vc
ẋ = 3P nvdisp −
rc
3P vdisp µa
p
= 3P nvdisp −
a
= 3P nvdisp − 3P nvdisp = 0 (177)
The ẏ component is given by:
~ c × ~rc ) + δ~rT (H
δ~v T (H ~˙ c × ~rc + H ~ c × ~vc )
ẏ =
||H~ c × ~rc ||
δ~rT (H~ c × ~rc )(H

~ c × ~rc )T (H ~˙ c × ~rc + H ~ c × ~vc )
−
||H~ c × ~rc ||3
T
~˙ c × ~rc + H

3P vdisp Ĥv × r̂c (H ~ c × ~vc )
T
= 3P nvdisp r̂c (Ĥc × r̂c ) −
~ c × ~rc ||
||H
T
~˙ c × ~rc + H

3P vdisp Ĥv × r̂c (Ĥc × r̂c )(H ~ c × ~rc )T (H ~ c × ~vc )
+
||H~ c × ~rc ||2
53
T
~˙ c × ~rc + H

3P vdisp Ĥv × r̂c (H ~ c × ~vc )
=−
||H~ c × ~rc ||
~ c × ~rc )T (H
3P vdisp (H ~˙ c × ~rc + H ~ c × ~vc )
+
||H~ c × ~rc ||2
T
~˙ c × ~rc + H

3P vdisp Ĥc × r̂c (H ~ c × ~vc )
=−
||H~ c × ~rc ||
3P vdisp (Ĥc × r̂c )T (H~˙ c × ~rc + H ~ c × ~vc )

+
||H~ c × ~rc ||
=0 (178)
Finally, the ż component can be found by:
~ c + δ~rT H
δ~v T H ~˙ c δ~rT H~ c (H
~ TH~˙ c )
c
ż = −
Hc Hc3
T T
~˙ c ~˙ c )

3P nvdisp r̂cT H~ c − 3P vdisp Ĥv × r̂c H 3P vdisp Ĥv × r̂c H ~ c (H
~ TH
c
= +
Hc Hc3
T T
~˙ c ~˙ c )

3P vdisp Ĥv × r̂c H 3P vdisp Ĥv × r̂c ~ c (H
H ~ cT H
= 3P nvdisp r̂cT Ĥc − + (179)
Hc Hc3
r̂c and Ĥc are orthogonal, so the first term goes to zero:
T T
~˙ c ~˙ c )

3P vdisp Ĥv × r̂c H 3P vdisp Ĥv × r̂c H ~ c (H
~ TH
c
ż = − +
Hc Hc3
T T
~˙ c ~˙ c )

3P vdisp Ĥv × r̂c H ~ cT H
3P vdisp Ĥc × r̂c Ĥc (H
=− + (180)
Hc Hc2
T
Ĥc × r̂c Ĥc = 0, so the 2nd term vanishes:
T
~˙ c

3P vdisp Ĥv × r̂c H
ż = −
Hc
54
T
3P vdisp Ĥv × r̂c ~rc × ~v˙ c
=−
Hc

3P vdisp (ĤvT ~rc )(r̂cT ~v˙ c ) − (ĤvT ~v˙ c )(r̂cT ~rc )
=− (181)
Hc
Ĥv = Ĥc and ~rc are orthogonal, therefore ĤvT ~rc = 0:
3P vdisp (ĤvT ~v˙ c )(r̂cT ~rc )

ż = (182)
Hc
~v˙ c = ~r¨c , and it is known from the two-body problem that the acceleration is in the
~ v and ~v˙ c are orthogonal. Therefore:
negative ~rc direction (Eq. (21)). In other words, H
ż = 0 (183)
In summary, The Hill frame state of the deputy, with respect to a chief satellite
deployed one orbit later, is:
 
0
 
−3P vdisp 
 
 
 
 0 
X̄dnew =
 
 (184)
 0 
 
 
 0 
 
0
So, the state of the deputy will be a y-offset with zero relative velocity (the deputy
will be in an ATO). There were a variety of assumptions that were employed in the
preceding development. First, the initial state of the deputy was propagated forward
one orbit using the HCW state transition matrix. This means that the result given
in Eq. (184) may not hold true in situations that violate the HCW assumptions
(described previously in the Background chapter). Two-body motion was assumed to
determine the angular velocity, meaning that the effects of nodal regression were not
55
considered. It was also assumed that the dispersal vehicle and the chief satellite lie
in the same orbit plane, but this assumption is reasonable as long as there is not a
significant cross-track component in the dispersal velocity. Additionally, there will be
small nonlinear effects that will cause the actual deputy state to deviate slightly from
Eq. (184). Note that there is no particular reason why the first deployed satellite
needs to be the deputy. If it was desirable for the second satellite to be the deputy,
the sign of Eq. (184) could simply be reversed (since the decision of a chief satellite
is arbitrary).
In order to quantify the effectiveness of the method, a simulation was run with
a dispersal vehicle in a 500 km altitude circular, equatorial orbit. The dispersal
velocity of the chief and deputy was assumed to be 1 m/s. The deputy initial state
was transformed into the ECI frame and propagated forward in time by one orbit
using the full nonlinear equations of motion (Eq. (186) without control). Then a
coordinate transformation was performed to determine the deputy state in a new
Hill frame centered on a chief satellite dispersed one orbit later. The result of the
simulation for the deputy state vector in the new Hill frame was:
 
−0.02 km
 
−17.07 km
 
 
 
 
 0 km 
X̄dnew =  (185)
 −0.014 m/s 
 
 
 −4
−3.03 × 10 m/s

 
0 m/s
The result predicted by Eq. (184) is y = −17.03 km, with 0 for the other Hill frame
components, so the differences between the predicted result and the simulation result
are relatively small in this case.
The method developed in this section shows that in the case of a circular chief
orbit and dispersals along the y-axis of the Hill frame, the relative velocity of the
56
satellites in a formation can be almost entirely eliminated by utilizing dispersals at
one orbit intervals. This approach has the advantage of reducing the amount of ∆V
required for initialization (since the satellites will not need to burn fuel to cancel
relative drift), while ensuring that collisions between satellites in the formation will
not occur.
3.6 Numerical Propagation of Solution
The trajectories are evaluated by propagating the orbits of the chief and deputy
satellite using absolute equations of motion including J2 perturbations and atmo-
spheric drag:
 
5( Zrdd )2 Xrdd

1−
µ 3 µ R 2 
e

~r¨d = − 3 ~rd − J2 2  1 − 5( Zr d )2 Yr d 
 
rd 2 rd rd  d

d
Zd 2 Zd
3 − 5( rd ) rd
 
Ẋ − ωE Yd
1 CD S d
q  d 
2 2 2
− ζγd (Ẋd − ωE Yd ) + (Ẏd − ωE Xd ) + Żd Ẏd − ωE Xd  + ~uECI (186)
 
2 md  
Żd
 
5( Zrcc )2 Xrcc

1−
µ 3 µ R 2 
e

~r¨c = − 3 ~rc − J2 2  1 − 5( Zrcc )2 Yrcc 
 
rc 2 rc rc  
Zc 2 Zc
3 − 5( rc ) rc
 
Ẋ − ωE Yc
1 CD Sc
q  c 
− ζγc (Ẋc − ωE Yc )2 + (Ẏc − ωE Xc )2 + Żc2 Ẏc − ωE Xc  (187)
 
2 mc  
Żc
where ~uECI is the control vector in the inertial frame. These equations of motion
provide a computationally efficient and relatively accurate method for evaluating the
effectiveness of the analytically determined thrust profiles. As can be seen from Eq.
57
(187), the chief satellite is uncontrolled. In the case of impulsive-thrust maneuvering,
the ~uECI term is dropped from Eq. (186), and copies of identical equations of motion
are used for the chief and deputy. In this case, the deputy satellite state vector is
updated after each coasting arc to include the ∆V (Eq. (41)).
In the continuous-thrust case, the control is not constrained to be constant.

This means that the mass flow rate (ṁ) is also not constant. The equation for the
specific impulse of a rocket engine with varying thrust is defined as [50]:
Rt
0
F (t)dt
Is = Rt (188)
gref 0 ṁ(t)dt
where gref is the reference acceleration (9.8066 m/s2 ). This equation can be differen-
tiated to yield:
F (t)
ṁ(t) = (189)
Is gref
with thrust magnitude F defined as F (t) = −m(t)U (t):
−m(t)U (t)
ṁ(t) = (190)
Is gref
Equation (190) is a differential equation that can be numerically integrated to track

the mass of the spacecraft throughout the maneuver.
In the case of impulsive thrust, the ∆V of a maneuver is solved for directly. How-
ever, the continuous-thrust case solves for the control acceleration. The continuous-
thrust ∆V can be found by simply integrating the control acceleration magnitude
U (t) over the transfer time tm . This integration can be performed analytically us-
ing the mathematical expression for U (t) or numerically using the values of ||~u(t)||
calculated throughout the maneuver.
Z tm
∆V = U (t)dt (191)
0
58
3.7 Algorithm Overview
The impulsive-thrust and continuous-thrust algorithms were both set up based

on equinoctial elements. This was necessitated by a desire to build a general algorithm
that can be used for formation flying in a variety of chief satellite orbits, including
orbits that are circular and equatorial.
The first step is to define the orbit of the chief satellite. Once this is done, the
initial state of the deputy can be defined. These initial states of the chief and deputy
are then propagated for a specified initial coast period before any maneuvering begins.
For the case of formation initialization, the deputy satellite is given an initial state
in the Hill frame, simulating dispersal from a container such as the CSD (the state
given by Eq. (184) was chosen for the formation initialization results shown in this
thesis, but the algorithm does not require this choice). For formation reconfiguration,
the initial conditions cannot be defined in such a straightforward manner, unless the
desired orbit has simple initial conditions (such as an ATO). For the cases of PCO and
GCO reconfiguration, the initial conditions are defined using a set of equations from
Table 2 and Eq. (89). Once the initial conditions of the deputy satellite are obtained,
the ECI state of the chief and deputy are propagated forward using identical copies
of the equations of motion (Eqs. (186 - 187), with the ~uECI term dropped).
Once the initial ECI states have been propagated for the specified coast period,
the deputy spacecraft begins the maneuvering phase. For the impulsive-thrust case,
the user defines a number of impulses (N ) and a set of parameters to define the desired
target orbit (such as ρ and αI for a desired target PCO). The solver then iterates until
it reaches a minimum-fuel solution that satisfies Eqs. (94) and (95). The solver then
returns a length N − 1 vector of the coast durations, a N × 3 matrix of ∆~v values (in
the Hill frame), and a wait duration (∆tw ). The states of the chief and deputy are
then propagated forward by ∆tw , and afterwards the first ∆~v is applied. The ∆~v is
translated into the ECI frame and applied to the deputy ECI state vector. The ECI
states of the chief and deputy are then propagated forward by the first coast duration,
59
and afterwards the second ∆~v is applied. This process is repeated until the end of
the maneuver phase.
For the continuous-thrust case, the user must define a transfer time and a desired
final state in the Hill frame. If the ~u required to complete the maneuver in the
specified time exceeds spacecraft performance parameters, the user is warned that the
spacecraft cannot complete the maneuver in the specified time. The control vector
~u is then rotated into the ECI frame and added to the equations of motion of the
deputy (Eq. (186)), and the chief and deputy states are propagated forward by the
transfer time.
Once the maneuvering phase is over, the states of the chief and deputy are
propagated forward by a specified time. If there is no formation maintenance, the
final ECI states of the chief and deputy are simply propagated using identical copies
of equations of motion with no control. In the case of impulsive-thrust formation
maintenance, the final states of the chief and deputy are propagated forward by the
DLQR time step size. The reference trajectory is then sampled, and the feedback
control law (Eq. (103)) is used to calculate the ∆~v . This ∆~v is then rotated into the
ECI frame and applied to the deputy state vector. The state of the chief and deputy
are then propagated forward by the time step, and the process is repeated until the
specified time for the final orbit propagation expires.
The continuous-thrust case for formation maintenance is virtually identical to

the case of continuous-thrust maneuvering. A transfer time (ttrans ) is specified, and
the reference trajectory is sampled at t = tcurrent + ttrans . The states of the chief
and deputy are then propagated forward by the transfer time using Eqs. (186-187).
Following the maneuver, the actual state of the deputy is updated (X̄factual and X̄f will
not be identical, since linearized dynamics are used in the thrust profile calculation),
and the process is repeated until the specified final orbit propagation time expires.
60
IV. Results
This chapter will present simulation results of the algorithms that have been outlined
in the previous chapter.
4.1 Formation Initialization
This section will present results for cases involving formation initialization. The
spacecraft coasts for a specified period simulating dispersal from a satellite dispenser
such as the CSD. Then the spacecraft executes the formation initialization maneuver.
The “+” marks the initial location of the deputy, and the “×” marks the location of
the chief (defined as the origin). The chief satellite was placed in a 500 km altitude
circular, equatorial orbit.
4.1.1 Impulsive-Thrust Initialization. For the impulsive-thrust formation

initialization results, the initialization geometry is shown on the left, and the right
hand side of the plots show the relative motion versus time. The circles mark the
locations of the impulses. A dashed vertical line marks the end of the coast phase,
and the dotted vertical lines mark the times of the impulses.
In Phase 1, the spacecraft executes an open-loop maneuver to the targeted

final state. In Phase 2, the orbit of the spacecraft is propagated with the formation
maintenance algorithm activated. This is to ensure that the spacecraft reaches the
desired final orbit, since linearization errors in the targeting dynamics can cause the
spacecraft to miss the desired final state at the end of the open-loop maneuver. The
∆V of Phase 1 and Phase 2 are summed to get the total initialization ∆V . This
approach is then compared to an alternate method where feedback is introduced after
each impulse.
Consider an impulsive-thrust formation initialization scenario involving a for-

mation of two 6U CubeSats with a mass of 10 kg. Assuming a 2-springs configuration,
the ejection velocity from the CSD is estimated to be 0.9 m/s (Fig. 8). The desired
formation is a ρ = 1 km PCO. The initialization maneuver with a one-orbit dispersal
61
interval (the deputy is dispersed one orbit after the chief), N = 2 and ∆tmax = 3
orbits is shown in Fig. 9. The ∆V for this maneuver was 1.94 m/s. It is clear from
the figure that linearization errors play a significant role in formation initialization,
as the > 14 km initial separation between chief and deputy is large enough to some-
what violate the d << r assumption of linearized dynamics. However, the open-loop
maneuver in Phase 1 gets the deputy sufficiently close for the closed-loop formation
maintenance algorithm to achieve the desired relative trajectory.
Figure 9: PCO Initialization Maneuver With N = 2, ∆tmax = 3 Orbits
Now it is desirable to see how varying maneuver parameters such as N and

∆tmax will affect the ∆V . Consider the same maneuver but with a more restrictive
∆tmax , shown in Fig. 10. The ∆V for this case was 2.92 m/s. Clearly, the allowable
transfer time has a significant effect on the ∆V for a formation initialization maneuver.
For the case in Fig. 11, the ∆tmax was restored to 3 orbits, and N was set to 3. In
this case the ∆V was 1.79 m/s. This shows that a marginal improvement in ∆V can
be made by increasing the number of burns. This is extrapolated even further in the
5-impulse, ∆tmax = 5 orbits case shown in Fig 12, with a ∆V of 1.73 m/s. Even
62
though the extended time for the open-loop maneuver causes it to be less accurate
than the previous cases, there is still an improvement in ∆V . This indicates that
there is likely a breaking point where the open-loop dynamics is too inaccurate, and
the total ∆V will increase due to the Phase 2 control having to correct the errors.
This is demonstrated in Fig. 13, where the ∆tmax was increased to 8 orbits, and the
∆V increased to 1.91 m/s. However, these errors can be corrected by introducing
feedback control in the Phase 1 maneuver. This case is shown in Fig. 14, where
the remaining impulses are re-computed after each burn. The ∆V was 1.31 m/s in
this case. Clearly, feedback control is necessary for the fuel-optimal maneuver in
most cases. Additionally, the Phase 2 curve in Fig 14 lies on top of the final orbit,
indicating that corrections to the Phase 1 maneuver were not necessary. This result
also indicates that higher-N maneuvers are likely to more accurately reach the desired
final orbit the first time, as there is more opportunity for feedback control.
Figure 10: PCO Initialization Maneuver With N = 2, ∆tmax = 1 Orbit
4.1.2 Low-Thrust Initialization. For the low-thrust initialization results,

the left side of the figures show the initialization geometry, and the right side shows
63
64
Figure 14: PCO Initialization Maneuver With N = 5, ∆tmax = 8 Orbits, Closed-

Loop Feedback
65
the control history for the maneuver. The horizontal bar indicates the control limit
for the spacecraft. This line is typically located at the top of the graph.
The low-thrust initialization maneuvers were performed with the spacecraft pa-
rameters defined in Table 3. These parameters were chosen to be representative of
a 6U CubeSat with an electrospray thruster. Like in the impulsive case, a one orbit
dispersal interval was implemented. It was found that feedback control is necessary
for most cases involving low-thrust formation initialization. For the results shown, the
spacecraft was given a position update once every 1/2 orbit in order to re-evaluate the
remainder of the trajectory. The spacecraft defined in Table 3 required 4.5 orbits to
complete the 1 km PCO initialization maneuver (Fig. 15). The ∆V for this maneuver
was 1.62 m/s, with a fuel cost of approximately 1.3 ×10−3 kg.
Table 3: Spacecraft Parameters

m F Is
10 kg 1.0 mN 1300 s
Figure 15: Low-Thrust PCO Initialization Maneuver, ∆t = 4.5 Orbits
Like in the impulsive case, it is desirable to see what effect lengthening the
transfer time has on the ∆V . The same PCO initialization maneuver with the transfer
time extended to 8 orbits is shown in Fig. 16. The ∆V for this maneuver was
66
1.52 m/s with a fuel cost of approximately 1.2 ×10−3 kg. This result is a marginal
improvement over the case shown in Fig. 15. The maneuver was executed once again
with a transfer time of 20 orbits, resulting in a ∆V of 1.49 m/s. These results indicate
that lengthening the transfer time does reduce the ∆V for low-thrust initialization,
but only at an incremental rate.
Figure 16: Low-Thrust PCO Initialization Maneuver, ∆t = 8 Orbits
Figure 17: Low-Thrust PCO Initialization Maneuver, ∆t = 20 Orbits
67
4.2 Formation Reconfiguration
This section will present results for cases involving formation reconfiguration.
PCO, GCO, and ATO reconfiguration maneuvers are examined. The chief satellite
was placed in a 500 km altitude circular, equatorial orbit.
4.2.1 Impulsive-Thrust Reconfiguration. The formation reconfiguration case

for impulsive-thrust is broken into four phases. For the first phase, referred to as the
“coast”, the initial conditions of the deputy satellite are propagated for one orbit to
show the initial orbit geometry. Then, the maneuver phase begins. The maneuver
phase is broken up into two parts – the “wait” period and the transfer. In the wait
period, the spacecraft coasts for an additional period, determined by the optimizer
(Eq. (93)). Then the spacecraft executes an open-loop N -impulse fuel-optimal ma-
neuver. Following the maneuver, the final orbit is propagated for one orbit to show
the final formation geometry.
Consider the case of a two-impulse GCO reconfiguration maneuver, shown in

Fig. 18. The maximum allowable transfer time for this maneuver was set to 3 orbits,
and the resulting ∆V was 1.11 m/s. Note that there is some fluctuation with the total
relative distance at the end of the maneuver. Even though GCO orbits have constant
chief-deputy spacing for all time in the HCW model, perturbations and nonlinearities
cause the real trajectory to deviate somewhat from this property. The fluctuations
also exist for the initial coast orbit, but they are harder to discern because they have
a smaller magnitude. These fluctuations can be mitigated somewhat by a formation
maintenance algorithm.
Even though the maximum transfer time for the result in Fig. 18 was set to 3
orbits, the fuel-optimal maneuver was completed in 1 orbit. This indicates that there
may not be any benefit to long transfer times in the case of 2-impulse reconfiguration.
Now, consider the same reconfiguration scenario, but with 3 impulses. This result is
shown in Fig. 19. The transfer time is extended in this case, but the resulting ∆V is
still 1.11 m/s. The same maneuver is again shown for 4 impulses and 5 impulses in
68
Figure 18: 1 km → 2 km GCO Reconfiguration with ∆tmax = 3 orbits, N = 2
Figs. 20 and 21, both with a ∆V of 1.11 m/s. The 5-impulse case was repeated, but
with the maximum allowable transfer time increased to 5 orbits (Fig. 22). The ∆V for
this maneuver was still 1.11 m/s. These results indicate that the GCO reconfiguration
∆V is approximately constant with the transfer time and number of impulses. A
similar result was observed for PCOs, summarized in Table 4 (plots are shown in
Figs. 43-47 in Appendix B).
Table 4: 1 km → 2 km PCO Reconfiguration Maneuver Results

N ∆tmax (orbits) ∆V m/s
2 3 1.24
3 3 1.24
4 3 1.24
5 3 1.25
5 5 1.24
69
Consider a ρ = 2 km → 1 km reconfiguration maneuver (Fig. 23), the reverse

of the maneuver shown in Fig. 18. The ∆V for this maneuver was 1.11 m/s. Note
that the ∆V for this case is equal to the 1.11 m/s ∆V for the 1 km → 2 km case.
The ρ = 1 km → 3 km case is shown in Fig. 24. The ∆V for this maneuver was
70
2.22 m/s. This is double the cost for the case shown in Fig. 18. These results
indicated a potential approximate linear relationship between |∆ρ| and ∆V for GCO
reconfiguration maneuvers for a given semimajor axis (where |∆ρ| = |ρf − ρ0 |) . This
prompted an investigation to see if general trends could be drawn between ∆ρ and
∆V . Simulations where run for a variety of orbit regimes and inclinations to see if
trends changed based on the reference orbit. The result is shown in Fig. 25. It was
found that the fuel-optimal ∆V , normalized by the chief orbit circular velocity Vcirc ,
∆ρ
is an approximate linear function of ac
for all orbit regimes that were analyzed.
The approximate linear trend is shown more explicitly in Fig. 26, where all of
the simulation results are plotted as one data set. As can be seen from the figure, all
of the simulation results lie on the linear curve. The slope of the linear curve implies
∆V ∆ρ
a 1:1 relation between Vcirc
and ac
. This suggests the approximate relationship for
GCO reconfiguration:
∆V ∆ρ
≈ (192)
Vcirc ac
√
∆ρ µ
∆V ≈ Vcirc ≈ 3/2 ∆ρ (193)
ac ac
71
Figure 25: GCO Formation Reconfiguration Results Using Impulsive-Thrust
Now consider an ATO reconfiguration maneuver, where the spacecraft is ma-

neuvering from a yd = 1 km ATO to a yd = 2 km ATO, shown in Fig. 27. The
∆V for the maneuver in Fig. 27 was 0.12 m/s. Note that for a similar change in
relative distance, the ATO reconfiguration maneuver is inexpensive in comparison to
the GCO (or PCO) reconfiguration maneuvers. This is likely due to the fact that
72
Figure 26: Approximate Linear Relationship For Impulsive GCO Reconfiguration
the ATO reconfiguration maneuver does not require any cross-track motion. Orbit
plane change maneuvers, which are oriented along the same direction as the angular
momentum vector (the Hill frame z-direction), are known to be costly from a ∆V
standpoint. The same reconfiguration maneuver was then repeated, but with N = 3.
The three-impulse ATO reconfiguration maneuver is shown in Fig. 28. The ∆V for
the three-impulse case in Fig. 28 is 0.06 m/s. Note that unlike the case of PCO and
GCO reconfiguration, there is a noticeable improvement in the ∆V when the number
of impulses is increased. This could potentially be explained by the fact that in the
ATO reconfiguration case, the center of relative motion is changing (an ATO can be
viewed as a 2:1 relative motion ellipse centered at ρy = yd with ρx = ρz = 0). This
means that unlike the case of PCO or GCO reconfiguration, the deputy satellite can
use a smaller ∆V and drift to the new location. In some ways, the ATO reconfigura-
tion maneuvers are more similar to the formation initialization results than the GCO
or PCO reconfiguration maneuvers.
73
The 4-impulse and 5-impulse reconfiguration cases are shown in Figs. 29 and
30. The ∆V for the 4-impulse and 5-impulse maneuvers was 0.05 m/s and 0.07 m/s,
respectively. However, when the ∆tmax for the 5-impulse case was extended to 5
orbits (Fig. 31), the ∆V reduced to 0.04 m/s. This indicates that in the case of
ATO reconfiguration, more impulses and a longer transfer time will generally result
in a lower ∆V . However, there are certain cases where having a higher number of
impulses is not necessarily more efficient.
Figure 27: 1 km → 2 km ATO Reconfiguration with ∆tmax = 3 orbits, N = 2
74
Figure 29: 1 km → 2 km ATO Reconfiguration r with ∆tmax = 3 orbits, N = 4
4.2.2 Low-Thrust Reconfiguration. The 1 km → 2 km GCO reconfiguration

maneuver was executed for the spacecraft in Table 3, and is shown in Fig. 32. The
maneuver was completed in four orbits and the ∆V was 1.31 m/s, corresponding to
a fuel cost of 1.0 × 10−3 kg. As expected, the ∆V cost for low-thrust reconfiguration
is somewhat higher than that of the impulsive-thrust case. The 1 km → 3 km GCO
reconfiguration maneuver using low-thrust is shown in Fig. 33. This maneuver was
completed in 7 orbits and required 2.62 m/s of ∆V , corresponding to a fuel cost of
75
2.1 × 10−3 kg. This indicates that as in the case of impulsive-thrust reconfiguration,
there is an approximate linear relationship between the ∆V and ∆ρ of a maneuver
for a given semimajor axis.
Figure 32: 1 km → 2 km Low-thrust GCO Reconfiguration Maneuver
It is now desirable to vary the spacecraft parameters to see what effect they
have on low-thrust reconfiguration scenarios. A large SmallSat (500 kg) with 23 mN
of thrust and an Isp of 2000 s would take 8 orbits to complete the 1 km → 2 km GCO
reconfiguration maneuver (Fig. 34), and would require 3.3 × 10−2 kg of propellant.
The large SmallSat took approximately twice as long to complete the same maneuver
76
Figure 33: 1 km → 3 km Low-Thrust GCO Reconfiguration Maneuver
as the 10 kg CubeSat. This indicates that a requirement for fast transfer times using
low-power devices would drive the design to a smaller satellite, since the maximum
achievable acceleration is limited by the upper bound on thrust for technologies such
as ion thrusters.
Figure 34: 1 km → 2 km Low-Thrust GCO Reconfiguration Maneuver with 500

kg, 23 mN Thrust Satellite
On the other hand, faster reconfiguration maneuvers on the order of one or

two orbits would necessitate high-power devices for all but the smallest of satellites.
For example, the 10 kg CubeSat with 4 mN of thrust could complete the 1 km →
2 km GCO reconfiguration maneuver in one orbit, and the 1 km → 3 km GCO
77
reconfiguration maneuver in two orbits (Figs. 35-36). This level of performance for
a CubeSat using low-thrust devices currently necessitates the use of a miniaturized
Hall thruster [41]. However, it is likely that the power requirements for a Hall-effect
thruster would be prohibitive for a spacecraft of this size. This indicates that a
requirement for fast transfer times would require either a low-thrust spacecraft with
a high power-to-mass ratio, or a design that utilizes chemical propulsion. Note that
the reconfiguration maneuvers shown here are relatively large (1 to 2 km change in
ρ for each maneuver). If the formations are smaller, the thrust requirements become
less prohibitive. For example, the 10 kg CubeSat with 1 mN of thrust could complete
a 0.5 km → 1 km GCO reconfiguration in 2 orbits (Fig. 37).
Figure 35: 1 km → 2 km Low-Thrust GCO Reconfiguration Maneuver with 10 kg,

4 mN Thrust Satellite
While the analysis for low-thrust reconfiguration has so far been focused on
general circular orbits, similar conclusions can be made for projected circular orbits
as well. The results of PCO reconfiguration cases are summarized in Table 5 (plots
are shown in Figs. 48-53 in Appendix B). Note how the transfer times and ∆V s are
somewhat larger in the case of PCOs, due to the additional cross-track magnitude.
78
Figure 36: 1 km → 3 km Low-Thrust GCO Reconfiguration Maneuver with 10 kg,
Figure 37: 0.5 km → 1 km Low-Thrust GCO Reconfiguration Maneuver
Table 5: PCO Reconfiguration Using Low Thrust

m F (mN) Isp (s) ∆ρ (km) ∆t (orbits) ∆V (m/s) ∆m (kg)
10 1 1300 1 4 1.49 1.2 ×10−3
10 1 1300 2 8 2.98 2.3 ×10−3
500 23 2000 1 9 1.49 3.8 ×10−2
10 4 1300 1 1.2 1.49 1.2 ×10−3
10 4 1300 2 2 2.99 2.3 ×10−3
10 1 1300 0.5 2 0.75 5.8 ×10−4
4.3 Drag Considerations
Up until the point, the presented results have assumed that the spacecraft are of
identical design flying at the same attitude. This represents the optimal case where
79
the effects of differential drag are minimized. This section will present results for
spacecraft flying at different attitudes (or of different designs) to examine the effects
of differential drag on satellite formations.
Consider the case shown in Fig. 14, where a 1 km PCO initialization maneuver
was executed resulting in a ∆V of 1.31 m/s. Now consider the same case, but this
time the deputy spacecraft has a projected area of 600 cm2 and the chief has a
projected area of 200 cm2 . This represents a worst-case differential drag scenario
for 6U CubeSats. A situation like this could arise if the chief spacecraft goes into
emergency mode while the deputy is maneuvering, or if the chief has some task to
perform that requires it to have a different attitude. The result using feedback control
is shown in Fig. 38, with a ∆V of 1.56 m/s. This value is a 19% increase in ∆V
from the ideal drag case. As can be seen from Fig. 38, the final PCO is also slightly
off-center (the middle of the 2:1 ellipse is not quite on the ×). This can be corrected
by increasing the number of impulses. The 6-impulse case is shown in Fig. 39, with
a ∆V of 1.39 m/s. As can be seen in Fig. 39, the accuracy of the final orbit can
be improved by increasing the number of impulses. Also, in this case, increasing the
number of impulses reduced the ∆V .
The low-thrust formation initialization cases shown previously were repeated as

worst-case differential drag scenarios, shown in Figs. 40-42. The ∆V s for the maneu-
vers in Figs. 40-42 were 1.62 m/s, 1.52 m/s, and 1.49 m/s, respectively. Interestingly,
these ∆V values are the same as those for the ideal drag cases shown in Figs. 15-17.
Clearly, differential drag does not have as significant of an effect on the low-thrust
algorithm, as compared to the impulsive-thrust case. This may be due to the fact that
the low-thrust initialization algorithm implements feedback on regular 1/2 orbit in-
tervals, whereas the impulsive-thrust algorithm often has long coast periods between
impulses. The long coast periods allow differential drag to build up a cumulative
effect that must be compensated for with the remaining burns of the maneuver.
80
Figure 38: 1 km PCO Initialization Maneuver with Worst-Case Differential Drag
Figure 39: 1 km PCO Initialization Maneuver with Worst-Case Differential Drag,

6 Impulses
81
using Low-Thrust, ∆t = 4.5 orbits

using Low-Thrust, ∆t = 8 orbits
82
using Low-Thrust, ∆t = 20 orbits
83
V. Conclusion
This chapter will present the contributions of this thesis and the key findings of the
simulation results in Chapter IV, following by a discussion of the limitations of the
algorithms. Finally, recommendations will be made for future research projects.
5.1 Contributions and Key Findings
This thesis presented the derivation of the GCO initial conditions in terms of
equinoctial elements, making it possible to implement GCO formations in equatorial
orbits. Physical significance of the bounded relative motion parameters is presented
for the case of general circular orbits. Impulsive-thrust and low-thrust algorithms
were developed for formation initialization, reconfiguration, and maintenance. The
1
algorithms are valid in circular and equatorial orbit regimes where singularities
are typically encountered. The trajectories are propagated with inertial equations
of motion, so the full nonlinear dynamics and J2 effects are captured, as well as
aerodynamic perturbations. Maneuver costs determined by the impulsive-thrust and
low-thrust algorithms were found to agree with similar analyses in the literature [4]
[15], giving a high degree of confidence to their accuracy. Methods were presented to
mitigate secular drift in the context of formation initialization. An approximate 1:1
∆V ∆ρ
relationship between Vcirc
and ac
was found for GCO reconfiguration using impulsive-
thrust.
It was found that open-loop control is adequate for many scenarios involving
formation reconfiguration. Results indicate that the formation reconfiguration cost
for PCOs and GCOs is approximately constant with the transfer time and number of
impulses. Low-thrust reconfiguration results show that requirements of fast transfer
times on the order of one or two orbits would drive the design of small satellites with
high power-to-mass ratios.
Formation initialization simulations demonstrate that longer transfer times and

a higher number of impulses generally reduces ∆V , but not in every case. Unlike
1
Code is available upon request. Contact email is Robert.LaRue@afit.edu
84
the case of formation reconfiguration, it was found that feedback control is necessary
for most cases involving formation initialization. An examination of worst-case drag
effects found that a significant increase in the impulsive-thrust ∆V can occur when
the spacecraft in a formation are flying at dissimilar attitudes for extended periods
of time. Differential drag effects on the low-thrust algorithm were less pronounced,
likely due to the regular interval of feedback in the low-thrust case.
5.2 Limitations of the Algorithms
One limitation of the developed algorithms is the use of the Schweighart-Sedwick

state transition matrix for maneuver targeting. While this STM accounts for J2
perturbations, it does not account for chief eccentricity. Therefore, the algorithms will
not be able to effectively generate trajectories for elliptical chief orbits. Additionally,
there are other relative motion models in the literature that are known to have higher
fidelity [27]. This indicates that it is likely that more optimal maneuvers could be
determined by using a higher-fidelity model. In the case of the impulsive thrust
algorithm, a state transition matrix that accommodates chief eccentricity could be
used, such as the Yamanaka-Ankerson [51] or Gim-Alfriend [14] models. For the
low-thrust algorithms, any model used would need to accomodate the assumptions
outlined in the Methodology. Cho and Park show explicitly in their article that a
variety of relative motion models (including Yamanaka-Ankerson) satisfy their method
[16].
Another notable limitation is the implementation of an exponential density

model in the calculation of the drag perturbations. This model does not take into ac-
count factors such as solar weather and earth geomagnetic activity, both of which the
atmospheric density is highly dependent on. While the implemented model provides
a reasonable approximation of differential drag effects, the accuracy of the estimation
could be significantly improved by implementing a higher-fidelity density model.
The algorithms account for J2 and drag perturbations, but do not consider other
perturbations such as solar radiation or third-body effects. This is sufficient for low-
85
Earth orbit, but these additional perturbations would need to be considered in order
for the algorithms to be accurate in other orbit regimes, such as the geostationary
belt. However, it is arguable that the spacecraft would have perturbations of a similar
magnitude.
An additional caveat is that perfect position knowledge and perfect control was
assumed for the work in this thesis. In reality a spacecraft’s position cannot be known
exactly, and thrusters have an associated error. These effects could be modeled by
adding an appropriate amount of random noise to spacecraft position and control in
the simulation.
5.3 Recommendations for Future Work
While this thesis took a step towards addressing the formation flying problem
for small satellites, there are many areas that could be improved or expanded upon.
As discussed in the previous section, the accuracy of the algorithms could be improved
by using a higher-fidelity relative motion model for targeting, and atmospheric density
estimation.
This thesis considered both impulsive and low-thrust methods, but only imple-
mented a single approach for each. The chosen methods were based somewhat on
the author’s experience and are not necessarily the most effective methods for Small-
Sat formation flying. There are a variety of formation control approaches that have
been derived for formation flying, many of which have been cited in the Introduction.
These could be adapted to similar problems as those addressed in this thesis. The
transfer times for the low-thrust algorithm were iterated manually to find trajectories
that satisfied the spacecraft constraints, but it would be possible to automate this
process. The fidelity of the algorithms could be improved by accounting for finite
burn times in the chemical propulsion case, rather than assuming impulsive burns.
The algorithms in this thesis compute fuel-optimal maneuvers, but there are other
cases (such as minimum time) that could be examined.
86
An analysis could be done to determine the “breaking point” of using relative
motion models for formation initialization targeting (i.e., how far away can the deputy
drift before linearization errors become so large that switching to a nonlinear targeting
approach is necessary?). The ∆V results for ATO reconfiguration were very low,
indicating that it is likely feasible to control and reconfigure these formations using
propellantless control methods such as solar sails or drag panels. The formation
initialization approach in this thesis assumed dispersal from the CSD, but there are a
variety of other SmallSat dispensers that could be considered. Finally, the algorithms
could be improved by accounting for attitude dynamics and position uncertainty.
87
Appendix A. State Transition Matrix
The Schweighart-Sedwick STM can be expressed as:
 
Φ11 Φ12 Φ13 Φ14 Φ15 Φ16
 
Φ21 Φ22 Φ23 Φ24 Φ25 Φ26 
 
 
 
Φ31 Φ32 Φ33 Φ34 Φ35 Φ36 
Φ=



Φ41 Φ42 Φ43 Φ44 Φ45 Φ46 
 
 
Φ51 Φ52 Φ53 Φ54 Φ55 Φ56 
 
Φ61 Φ62 Φ63 Φ64 Φ65 Φ66
The nonzero entries of the matrix are:
1
Φ11 = (4(1 + s) − (3 − 5s) cos (ξt))
1−s
1
Φ14 = √ sin (ξt)
n 1−s
√
2 1+s
Φ15 = (1 − cos (ξt))
n(1 + s)
√
2 1 + s(3 + 5s)
Φ21 = √ (sin (ξt) − ξt)
(1 − s) 1 − s
Φ22 = 1
√
2 1+s
Φ24 = (cos (ξt) − 1)
n(1 − s)

1 4(1 + s)
Φ25 = √ sin (ξt) − (3 + 5s)t
1−s n 1−s
Φ33 = cos (qt)
1
Φ36 = sin (qt)
q
n(3 + 5s)
Φ41 = √ sin (ξt)
1−s
Φ44 = cos (ξt)
88
√
2 1+s
Φ45 = √ sin (ξt)
1−s
√
2n 1 + s(3 + 5s)
Φ51 = (cos (ξt) − 1)
1−s
√
2 1+s
Φ54 = − √ sin (ξt)
1−s
1
Φ55 = (4(1 + s) cos (ξt) − (3 + 5s))
1−s
Φ63 = −q sin (qt)
Φ66 = cos (qt)
The parameters are:
3J2 Re2
s= (1 + 3 cos (2ic ))
8rc2
√
ξ =n 1−s
√
c= 1+s
r
µ
n=
rc3
δh0 = hd − hc

cot (ic ) sin (id ) − cos (id ) cos (δh0 )
γ0 = acot
sin (δh0 )
3 J2 nRe2
ḣd = − cos (id )
2 rc2
3 J2 nRe2
ḣc = − cos (ic )
2 rc2
q = nc − (cos (γ0 ) sin (γ0 ) cos (δh) − sin (γ0 )2 cos (id ))(ḣd − ḣc ) − ḣd cos id
if id ≈ ic
3J2 nRe2
q = nc + 2
cos (ic )2
2rc
89
Appendix B. Additional Results
Figure 43: 1 km → 2 km PCO Reconfiguration with ∆tmax = 3 orbits, N = 2
90
91
Figure 48: 1 km → 2 km Low-Thrust PCO Reconfiguration Maneuver
92
Figure 49: 1 km → 3 km Low-Thrust PCO Reconfiguration Maneuver
Figure 50: 1 km → 2 km Low-Thrust PCO Reconfiguration Maneuver with 500 kg,

93

94
Figure 53: 0.5 km → 1 km Low-Thrust PCO Reconfiguration Maneuver
95
Bibliography
1. Heidt, H., Puig-Suari, J., Moore, A. S., Nakasuka, S., and Twiggs, R., “CubeSat:
A New generation of picosatellite for education and industry low-cost space ex-
perimentation,” 14th Annual AIAA/USU Conference on Small Satellites, Logan,
UT, 2000.
2. Bouwmeester, J. and Guo, J., “Survey of worldwide pico- and nanosatellite mis-
sions, distributions and subsystem technology,” Acta Astronautica, Vol. 67, 2010,
pp. 854 – 862.
3. Gill, E., Sundaramoorthy, P., Bouwmeester, J., Zandbergen, B., and Reinhard,
R., “Formation flying within a constellation of nano-satellites: The QB50 mis-
sion,” Acta Astronautica, Vol. 82, 2013, pp. 110 – 117.
4. Alfriend, K. T., Vadali, S. R., Gurfil, P., How, J. P., and Breger, L. S., Spacecraft
Formation Flying, Elsevier Astrodynamics Series, Elsevier, Oxford, UK, 2010.
5. Tooley, C. R., Black, R. K., Robertson, B. P., Stone, J. M., Pope, S. E., and
Davis, G. T., “The Magnetospheric Multiscale Constellation,” Springer Science,
Vol. 199, 2016, pp. 23–76.
6. Queen, S., Shah, N., Benegalrao, S., and Blackman, K., “Generalized Momentum
Control of the Spin-Stabilized Magnetospheric Multiscale Formation,” AAS/A-
IAA Astrodynamics Specialist Conference, Vol. 15, Vail, CO, 2015, pp. 1–17.
7. D’Amico, S., Ardaens, J. S., and DeFlorio, S., “Autonomous Formation Flying
Based on GPS - PRISMA Flight Results,” Acta Astronautica, Vol. 82, 2013,
pp. 69– 79.
8. Chung, J. S., Bandyopadhyay, S., Foust, R., and Subramanian, G. P., “Review of
Formation Flying and Constellation Missions Using Nanosatellites,” Journal of
Spacecraft and Rockets, Vol. 53, No. 3, 2013.
9. Roscoe, C., Westphal, J., Lutz, S., and Bennet, T., “Guidance, Navigation, and
Control Algorithms for CubeSat Formation Flying,” 38th AAS Guidance and
Control Conference, Vol. 15, Breckenridge, CO, 2015.
10. Bonin, G., Roth, N., Armitage, S., Newman, J., Risi, B., and Zee, R. E., “CanX-4
and CanX-5 Precision Formation Flight: Mission Accomplished!” 29th Annual
AIAA/USU Conference on Small Satellites, Logan, UT, 2015, pp. 1–15.
11. Vaddi, S. S., Alfriend, K. T., Vadali, S. R., and Sengupta, P., “Formation Es-
tablishment and Reconfiguration Using Impulsive Control,” Journal of Guidance,
Control, and Dynamics, Vol. 28, No. 2, 2005, pp. 262–268.
12. Palmer, P., “Optimal Relocation of Satellites Flying in Near-Circular Orbit For-
mations,” Journal of Guidance, Control, and Dynamics, Vol. 29, No. 3, 2006,
pp. 519–526.
96
13. Vignal, P. and Pernicka, H., “Low-Thrust Spacecraft Formation Keeping,” Jour-
nal of Spacecraft and Rockets, Vol. 43, No. 2, 2006, pp. 466–475.
14. Gim, D. and Alfriend, K. T., “State Transition Matrix of Relative Motion for
the Perturbed Nonsingular Reference Orbit,” Journal of Guidance, Control, and
Dynamics, Vol. 26, No. 6, 2003, pp. 956–971.
15. Yan, H. and Alfriend, K. T., “Approximate Minimum Energy Control Laws for
Low-Thrust Formation Reconfiguration,” Journal of Guidance, Control, and Dy-
namics, Vol. 30, No. 4, 2007, pp. 1182–1185.
16. Cho, H. C. and Park, S. Y., “Analytic Solution for Fuel-Optimal Reconfiguration
in Relative Motion,” Journal of Optimization Theory and Applications, Vol. 141,
No. 3, 2009, pp. 495–512.
17. Lee, S. and Park, S. Y., “Approximate Analytical Solutions to Optimal Reconfig-
uration Problems in Perturbed Satellite Relative Motion,” Journal of Guidance,
18. Massari, M. and Bernelli-Zazzera, F., “Optimization of Low-Thrust Reconfig-
uration Maneuvers for Spacecraft Flying in Formation,” Journal of Guidance,
19. Cho, H. C., Park, S. Y., Yoo, S. M., and Choi, K. H., “Analytical Solution to
Optimal Relocation of Satellite Formation Flying in Arbitrary Elliptic Orbits,”
Aerospace Science and Technology, Vol. 25, No. 5, 2012, pp. 161–176.
20. Tschauner, J. and Hempel, P., “Rendezvous zu einem in elliptischer Bahn um-
laufenden Ziel,” Astronautica Acta, Vol. 11, No. 2, 1965, pp. 104–109.
21. “General Payload User’s Guide,” Spaceflight, Inc., 2015, SF-2100-PUG-00001
Rev. F.
22. “Canisterized Satellite Dispenser Data Sheet,” Planetary Systems Corporation,
2017, 2002337 Rev. E.
23. Hill, G. W., “Researches in the Lunar Theory,” American Journal of Mathemat-
ics, Vol. 1, No. 1, 1878, pp. 5–26, 129–147, 245–260.
24. Clohessy, W. H. and Wiltshire, R. S., “Terminal Guidance System for Satellite
Rendezvous,” Journal of the Aerospace Sciences, Vol. 27, 1960, pp. 653–658, 674.
25. Vallado, D. A., Fundamentals of Astrodynamics and Applications, Microcosm
Press, Hawthorne, CA, 4th ed., 2013.
26. Schweighart, S. A. and Sedwick, R. J., “High-Fidelity Linearized J2 Model
for Satellite Formation Flight,” Journal of Guidance, Control, and Dynamics,
Vol. 25, No. 6, 2002, pp. 1073–1080.
27. Sullivan, J., Grimberg, S., and D’Amico, S., “Comprehensive Survey and As-
sessment of Spacecraft Relative Motion Dynamics Models,” Journal of Guidance,
Control, and Dynamics, 2017, pp. 1–21.
97
28. Johnson, K. W., Approaches for Modeling Satellite Relative Motion, PhD disser-
tation, Texas A & M University, 2016.
29. Vadali, S. R., Sengupta, P., Yan, H., and Alfriend, K. T., “Fundamental Frequen-
cies of Satellite Relative Motion and Control of Formations,” Journal of Guidance,
30. Vadali, S. R., Vaddi, S. S., and Alfriend, K. T., “An Intelligent Control Concept
for Formation Flying Satellite Constellations,” International Journal of Robust
and Nonlinear Control , Vol. 12, 2002, pp. 97–115.
31. Montenbruck, O. and Gill, E., Satellite Orbits: Models, Methods, and Applica-
tions, Springer, New York, NY, 3rd ed., 2005.
32. Cefola, P. J., “Equinoctial Orbit Elements - Application To Artificial Satellite
Orbits,” AIAA/AAS Astrodynamics Conference, No. 72-937, 1972, p. 3.
33. Wiesel, W., Modern Astrodynamics, Aphelion Press, Beavercreek, OH, 2010.
34. Hintz, G. R., “Survey of Orbit Element Sets,” Journal of Guidance, Control, and
Dynamics, Vol. 31, No. 3, 2008, pp. 785–789.
35. Schaub, H. and Junkins, J. L., Analytical Mechanics of Space Systems, American
Institute of Aeronautics and Astronautics, Reston, VA, 3rd ed., 2014.
36. Markley, F. L. and Crassidis, J. L., Fundamentals of Spacecraft Attitude Deter-
mination and Control , Space Technology Library, Springer, New York, 2014.
37. Brouwer, D., “Solution of the Problem of Artificial Satellite Theory without
Drag,” The Astronomical Journal , Vol. 64, No. 1274, 1959, pp. 378–397.
38. Greenwood, D. T., Classical Dynamics, Dover Science Books, Dover, Mineola,
NY, 1997.
39. Lanczos, C., The Variational Principles of Mechanics, Dover Books on Physics,
Dover, Mineola, NY, 1986.
40. Gim, D. and Alfriend, K. T., “Satellite Relative Motion Using Differential
Equinoctial Elements,” Celestial Mechanics and Dynamical Astronomy, Vol. 92,
2005, pp. 295–336.
41. “Small Spacecraft Technology State of the Art,” Tech. Rep. TP2015216648,
NASA Ames Research Center, 2015.
42. Conway, B. A., Spacecraft Trajectory Optimization, Cambridge Aerospace Series,
Cambridge University Press, Cambridge, 2014.
43. Prussing, J. E. and Conway, B. A., Orbital Mechanics, Oxford University Press,
Oxford, 2013.
44. Wiesel, W. E. and Alfano, S., “Optimal Many-Revolution Orbit Transfer,” Jour-
nal of Guidance, Control, and Dynamics, Vol. 8, 1985, pp. 155–157.
98
45. Goebel, D. M. and Katz, I., “Fundamentals of Electric Propulsion: Ion and Hall
Effect Thrusters,” Jet Propulsion Laboratory Space Science and Technology Se-
ries, 2008.
46. Schaub, H. and Alfriend, K. T., “J2 Invariant Relative Orbits for Spacecraft For-
mations,” Celestial Mechanics and Dynamical Astronomy, Vol. 79, No. 2, 2001,
pp. 77–95.
47. Vaddi, S. S., Vadali, S. R., and Alfriend, K. T., “Formation Flying: Accommodat-
ing Nonlinearity and Eccentricity Perturbations,” Journal of Guidance, Control,
and Dynamics, Vol. 26, No. 2, 2003, pp. 214–223.
48. Vadali, S. R., Schaub, H., and Alfriend, K. T., “Initial Conditions and Fuel-
Optimal Control for Formation Flying of Satellites,” AIAA GNC Conference,
No. AIAA 99-4265, Portland, OR, 1999.
49. Kunz, D. L., Intermediate Dynamics for Aeronautics and Astronautics, Head-
master Press, Centerville, OH, 2015.
50. Sutton, G. P. and Biblarz, O., Rocket Propulsion Elements, John Wiley and Sons,
Hoboken, NJ, 8th ed., 2010.
51. Yamanaka, K. and Ankerson, F., “New State Transition Matrix for Relative Mo-
tion on an Arbitrary Elliptical Orbit,” Journal of Guidance, Control, and Dy-
namics, Vol. 25, No. 1, 2002, pp. 60–66.
99
Form Approved
REPORT DOCUMENTATION PAGE OMB No. 0704-0188
The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of
information, including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188),
1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any
penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.
PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.
1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To)
22-03-2018 Master's Thesis September 2016 - March 2018
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Algorithms for Small Satellite Formation Flying
5b. GRANT NUMBER
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER

LaRue, Robert B, 2LT 18YD334D
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

REPORT NUMBER
Graduate School of Engineering and Management (AFIT/EN) AFIT-ENY-MS-18-M-273
2950 Hobson Way
Wright-Patterson AFB OH 45433-7765
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S)
Air Force Research Laboratory AFRL/RV
Charlene Jacka , Deputy Program Manager
20170 E Kirtland AFB, 11. SPONSOR/MONITOR'S REPORT
Kirtland AFB, NM 87117 NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT

DISTRIBUTION A. APPROVED FOR PUBLIC RELEASE; DISTRIBUTION UNLIMITED
13. SUPPLEMENTARY NOTES

This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States.
14. ABSTRACT
This thesis presents algorithms for spacecraft formation flying using impulsive-thrust and low-thrust methods. The general circular
orbit formation initial conditions are derived in terms of equinoctial elements. Physical significance of the bounded relative motion
parameters is presented for the case of general circular orbits. The developed algorithms are posed in terms of equinoctial elements
for a singularity-free approach. The algorithms are assessed by numerical propagation of the inertial equations of motion with J2
and drag perturbations. Methods are presented for minimizing the ∆V required for formation initialization. An examination of the
performance of open-loop and closed-loop control is provided for formation initialization and reconfiguration. The effects of
differential drag on small satellite formations is analyzed. The developed algorithms are used to examine the trade space and
quantify how spacecraft design parameters affect formation flying scenarios.
15. SUBJECT TERMS

Astrodynamics, Relative Motion, Formation Flying, Formation Design, Guidance and Control
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF RESPONSIBLE PERSON
a. REPORT b. ABSTRACT c. THIS PAGE ABSTRACT OF Lt Col Kirk W. Johnson, AFIT/ENY
PAGES
UU 19b. TELEPHONE NUMBER (Include area code)
U U U 112
937-255-3636 x4285 kirk.johnson@afit.edu
Standard Form 298 (Rev. 8/98)
Prescribed by ANSI Std. Z39.18

Algorithms For Small Satellite Formation Flying

Uploaded by

Copyright:

Available Formats

Algorithms For Small Satellite Formation Flying

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Algorithms For Small Satellite Formation Flying

Uploaded by

Copyright:

Available Formats

Air Force Institute of Technology

Theses and Dissertations Student Graduate Works

Algorithms for Small Satellite Formation Flying

Follow this and additional works at: https://scholar.afit.edu/etd

Robert B. LaRue, 2nd Lieutenant, USAF

DEPARTMENT OF THE AIR FORCE

AIR FORCE INSTITUTE OF TECHNOLOGY

DISTRIBUTION A. APPROVED FOR PUBLIC RELEASE; DISTRIBUTION IS UNLIMITED.

Algorithms for Small

Presented to the Faculty

Robert B. LaRue, B.S.

DISTRIBUTION A. APPROVED FOR PUBLIC RELEASE; DISTRIBUTION IS UNLIMITED.

ALGORITHMS FOR SMALL

Robert B. LaRue, B.S.

Lt Col Kirk W. Johnson, PhD

William E. Wiesel, PhD

Capt Joshuah A. Hess, PhD

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

Appendix A. State Transition Matrix . . . . . . . . . . . . . . . . . . 88

Appendix B. Additional Results . . . . . . . . . . . . . . . . . . . . . 90

1. Exponential Density Model . . . . . . . . . . . . . . . . . . . . 17

A satellite formation is defined by NASA’s Goddard Space Flight Center (GSFC)

While there are benefits to using a formation of satellites as opposed to a single

There have been numerous studies on the reconfiguration of satellite relative

To summarize, the goals of the research presented are as follows:

2. Compare and evaluate a variety of formation types

3. Examine design tradespace to see how spacecraft parameters affect formation

This Introduction chapter will be followed by a Background chapter, which sum-

2.1 Relative Motion and Formation Flying

One of original models for satellite relative motion is the Hill/Clohessy-Wiltshire

2ẏ0 ẋ0 2ẏ0

2ẋ0 4ẏ0 2ẋ0

ẏ0 = −2nx0 (7)

2.2 Satellite Formations

x(t) = ρx sin (nt + αx ) (8)

y(t) = ρy + 2ρx cos (nt + αx ) (9)

z(t) = ρz sin (nt + αz ) (10)

αz = atan2(nz0 , ż0 ) (15)

Figure 2: PCO Relative Orbit Geometry at t = 0

2.2.2 Along-Track Orbits. Along-track orbits (ATOs) are a type of space-

2.3 Orbit Element Sets

where a is the semimajor axis, e is the eccentricity, i is the inclination, h is the

where λ = l + g, q1 = e cos g, and q2 = e sin g. It is important to note that the

2.4 Orbit Element Differences

2.5 Orbital Perturbations

The motion of a spacecraft is primarily governed by the gravity of the primary

2.5.1 Perturbations due to an Aspherical Gravity Field. The analysis pre-

2.5.2 Perturbations due to Atmospheric Drag. The perturbation due to

~vrel = ~v − (~ωE × ~r) (25)

where ωE = 7.2921158553 × 10−5 rad/s [36].

Table 1: Exponential Density Model

can be approximated as γ = r − RE . Then, Eq. (24) can be used to estimate the

2.6 Mean Orbit Elements

D̄ = (l, g, h, L, GH) (28)

where H1 is the perturbing Hamiltonian and  is the perturbation. In this case,

Figure 5: Categories of Small Satellites

Satellites of the sizes illustrated in Fig. 5 tend to lend themselves towards

2.8 State Transition Matrix

X̄t = Φt X̄0 (39)

where H1 is the perturbing Hamiltonian and is the perturbation. In this case,