Simultaneous Observation Association and Maneuver Reconstruction for Non-Keplerian

Initial Orbit Determination Using Nonlinear Programming

Casey R. Heidrich and Marcus J. Holzinger

University of Colorado Boulder, Boulder, CO USA

ABSTRACT

Accurate initial orbit determination from optical or range measurements presents numerous challenges for non-Keplerian systems.
In addition to inherent chaos and low predictability of multi-body orbits over long time scales, a spacecraft may also perform
unknown maneuvers between observations. Solving the measurement-to-measurement association problem with dynamic mis-
modeling becomes intractable without a priori knowledge of the target orbit. This work seeks to overcome these challenges by
leveraging implicit integration methods and nonlinear programming to perform initial orbit determination, maneuver reconstruc-
tion, and measurement association of active spacecraft. The approach greatly enhances operational capability for detection and
tracking of objects in cislunar space.

1. INTRODUCTION

Reliable initial orbit determination (IOD) is a critical capability for space situational awareness (SSA). The funda-
mental objective of IOD is the estimation of a full spacecraft state using electro-optical (EO) or range measurements.
Historical IOD methods are built upon analytical solutions from two-body Keplerian mechanics, such as Gauss’s
method or double-r [1]. These assumptions fail in complex multi-body systems where traditional two-body me-
chanics are invalid. More recent work has developed IOD methods for non-Keplerian systems using optimization
techniques [2, 3, 4]. These methods approach the limitations of classical IOD methods by approximating solutions of
nonlinear system trajectories using numerical integration. However, chaotic dynamics such as the circular restricted
three-body problem (CR3BP) become difficult to predict over long time scales, and these methods quickly become
unstable with large measurement gaps or insufficient observations.
Current research has focused on improvements to both the accuracy and reliability of IOD methods in complex orbit
environments [5]. The algorithm developed in [6] demonstrates a novel approach to non-Keplerian IOD using collo-
cation methods. The methodology involves transcription of a continuous-time optimization problem to a large-scale
sparse nonlinear programming (NLP) problem. These methods are considered implicit integrators that enforce sys-
tem dynamics through collocation constraints (in contrast to time-marching algorithms such as explicit Runge-Kutta
schemes [7]). The collocation IOD approach provides many benefits, including a large region of convergence from
poor initial guesses, as well as improved stability over long observing gaps. Whereas prior work focused on quiescent
(non-maneuvering) objects, operational realities will require greater flexibility in terms of dynamical modeling and
assumptions for IOD of maneuvering spacecraft.
The measurement association problem is an open area of research in SSA, the purpose of which is to identify obser-
vations corresponding to the same object at different epochs. The challenge of associating uncorrelated observations
has been approached using covariance-based methods [8], Bayesian inference [9], and admissible region (AR) based
optimization methods [10], among others. The latter involves the solution of a constrained optimization problem for
undetermined quantities satisfying AR constraints at each measurement. Similarly, maneuver reconstruction is an im-
portant topic in SSA for inferring spacecraft behavior from sparse measurement information. Classical sequential or
batch estimation algorithms often break down when the true system behavior does not match an assumed dynamics
model. Prior works use control distance metrics for maneuver detection and object correlation by solving an uncertain
two-point boundary value problem (BVP) for the optimal control policy between observations [11, 12].
Approved for public release; distribution is unlimited. Public Affairs release approval AFRL-2024-4378. The views expressed are those of the
authors and do not reflect the official guidance or position of the United States Government, the Department of Defense or of the United States Air
Force.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
This work develops a novel non-Keplerian IOD algorithm incorporating maneuver reconstruction and observation
association without prior information. Problem objectives and constraint structures are considered. The overall prob-
lem structure is retained such that efficient NLP solvers may take advantage of inherent constraint sparsity for rapid
solution. Results demonstrate capability for reconstructing large maneuvers of cislunar spacecraft with limited infor-
mation, such as angles-only measurements across observation gaps of days or weeks. The algorithm produces accurate
results without a close initial guess for the spacecraft state or control history. For example, convergence is shown for
initialization at one of the Earth-Moon Lagrange points and null (all zeros) for the control history. This work also
develops metrics to quantify the accuracy and consistency of IOD solutions.
Robust and adaptable IOD methods are critical to continued success of operations in cislunar space. While the as-
sumptions of this work do not preclude applications in traditional two-body dynamics (such as near-Earth orbits),
they are general enough to allow for rapid IOD with complex multi-body systems and maneuvering targets. Further-
more, including admissible region information as NLP constraints provides a direct approach for implicit measurement
association. Methods for solution with representative applications are presented to validate the algorithm.

2. BACKGROUND

The following section describes the relevant background material and dynamical models for the IOD algorithm de-
veloped in this paper. The spacecraft state is described by its position and velocity coordinates in a rotating frame
coincident with the synodic period of the Earth-Moon system
 T
x = x y z ẋ ẏ ż (1)

The problem dynamics are described by the circular restricted three-body problem (CR3BP) [13]

ẍ = 2ẏ + x − (1 − µ)(x + µ)r1−3 − µ(x − 1 + µ)r2−3 + u1 (2a)

ÿ = −2ẋ + y − (1 − µ)yr1−3 − µyr2−3 + u2 (2b)
z̈ = −(1 − µ)zr1−3 − µzr2−3 + u3 (2c)

where u1 , u2 , and u3 represent external control inputs (thrust) to the system. The r1 and r2 quantities represent distances
of the spacecraft relative to the primary and secondary bodies, respectively. In the absence of external control inputs,
this system gives rise to the pseudo-potential function
1 1−µ µ
U = (x2 + y2 ) + + (3)
2 r1 r2
which in turn admits the well-known Jacobi constant

C = 2U − v2 (4)

where
v2 = ẋ2 + ẏ2 + ż2 (5)

Although the CR3BP dynamics do not consider more realistic four-body effects or other perturbations, these equations
are generally considered sufficient for describing the motion of cislunar orbits for mission design purposes.
A generic observing scenario in the CR3BP coordinate system is illustrated in Fig. 1. Angles measurements are
parameterized in terms of the azimuth α and declination δ angles of the relative position vector. If we define r ∈ Rn/2
and o ∈ Rn/2 as the target and observer positions, respectively, then the relative position vector is l = r − o . It is
straightforward to define the angles measurements as
 
−1 ly
α = tan (6a)
lx
 
lz
δ = sin−1  q  (6b)
2 2 2
lx + ly + lz

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 Fig. 1: Depiction of instantaneous angle measurements in rotating frame coordinates (not to scale).

However, to avoid discontinuities with the arctangent operation it is often desirable to use a line-of-sight (LOS) vector
measurement model [14] defined as  
cos α cos δ
l̂l =  sin α cos δ  (7)
sin δ
With angles and angle rates measurements generated from optical tracklets, the position and velocity states of the
target can be inferred as
r = o + ρ l̂l (8)
and
ṙr = ȯo + ρ̇ l̂l + ρ α̇ l̂l α + ρ δ̇ l̂l δ (9)
where ρ is the target range, ρ̇ is the range rate, and
   
− sin α cos δ − cos α sin δ
l̂l α =  cos α cos δ  , l̂l δ =  − sin α sin δ  (10)
0 cos δ

The problem of measurement association is common when one or more uncorrelated tracks (UCT) cannot be cor-
related to any known object. This problem has been studied extensively using covariance-based track association
methods [8, 15] and admissible regions [16, 10, 17]. This paper focuses on the latter, as covariance-based methods
require both a full state and uncertainty estimate, which are not available prior to successful IOD. For two or more
UCT observations to be correlated, the admissible region approach requires range and range-rate values to satisfy
orbit constraints. With Keplerian systems, it is common to define an admissible region based on orbital energy and
eccentricity. However, these quantities are not applicable to non-Keplerian dynamics. Measurement association for
cislunar orbits can leverage analogous constraints on the Jacobi constant U and two-body orbital energy ε of the target
relative to the primary and secondary bodies [18, 19].
An example of admissible region constraints for a cislunar optical observing scenario is shown in Fig. 2. The intersec-
tion of the interior region of these constraints defines the possible combinations of range and range-rate at the specified
measurement epoch. Measurement association using admissible regions typically relies on forward or backward nu-
merical integration of a hypothesized range and range-rate pair to another measurement. This task becomes arduous
with an increasing number of measurements. In the following sections, measurement association is achieved using
admissible region methodologies by implicitly enforcing state constraints at each measurement epoch.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
Fig. 2: Example of admissible region constraints in range and range-rate coordinates for an L1 Northern Halo observer.
Jacobi constant bounds are based on limits of the L2 Southern Halo family.

3. METHODOLOGY

The contribution of this work is a collocation-based IOD framework for non-Keplerian systems. The novel approach
to this problem leverages implicit integration of the system dynamics, with relevant constraints included in an NLP
framework. The approach is highly robust and can converge to accurate solutions without a close initial guess for the
initial state and control inputs. First, the time domain is discretized into nodes (or “knots”), which serve as collocation
points in the solver. A transcription method is applied to convert the continuous-time control problem into a parameter
optimization problem. Control inputs are designed to minimize the overall control effort and fuel necessary to achieve
the constructed orbit solution. The method is also modified to enable IOD with impulsive maneuvers, which are more
typical of spacecraft station-keeping. The following section outlines this process in more detail.
3.1 Problem Formulation
Motion of the system in (2) is described by states x ∈ Rn and control inputs u ∈ U , where U is a compact and bounded
set in Rm .
ẋx(t) = f (xx(t), u (t),t) (11)
The measurement model of an observation y ∈ Rq at time ti is modeled as

y i = h (xx(ti ),ti ) + v i (12)

where v ∈ Rq is a zero-mean Gaussian white noise process with intensity E[vv(ti )vv(ti )T ] = R i ∈ Sq++ . Suppose a
collection of optical measurements {yy1 , y 2 , . . . , y p } are collected by an observer. The IOD objective is to reconstruct
the state trajectory and possible control inputs producing the observations at each measurement epoch, subject to
admissible region constraints. In [6], this problem is summarized by a weighted least squares cost
p
J = ∑ (yyi − h (xx(ti ),ti ))T R −1 yi − h (xx(ti ),ti ))
i (y (13)
i+1

The solution to this objective is known to provide the maximum likelihood estimate (MLE) under certain regularity
assumptions for linear systems [20]. In the sections to follow, we will relax this objective by defining the measurement
residuals as
ζ i = y i − h (xx(ti ),ti ) (14)

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
and optimizing for control effort instead. Let us also define an augmented vector of the measurement residuals and
admissible region constraints at each measurement epoch
 T
c (xx(ti ),ti ) = ζ i Ci ε1,i ε2,i (15)

where ε1 and ε2 are the two-body orbital energies of the spacecraft with respect to the primary and secondary bodies,
respectively. Next, we consider control distance metrics for modeling spacecraft maneuvers [9, 21]. Track association
methods using control distances model the likelihood of a maneuver between two UCT by minimizing the performance
cost Z tf
P= u (τ)T u (τ) dτ (16)
t0

as the distance between an initial state (xx0 ,t0 ) and final state (xx f ,t f ). This energy-like objective minimizes the integral
of sum-squared control effort over the trajectory, which is typical of low-thrust continuous maneuvers.
The complete IOD problem structure is outlined in (17), which includes an integrated control effort cost and observa-
tion constraints.  R tf T
 minimize t0 u (τ) u (τ) dτ
x ,u ∈U
u



Problem P(t) subject to: ẋx(t) − f (xx(t), u (t),t) = 0 (17)


 c L ≤ c (xx(ti ),ti ) ≤ cU
for i = 1, . . . , p


The quantities cL , and cU represent bounds on the measurement constraints in (15). In stated form, this problem is
challenging to solve with continuous-time dynamics and point constraints at each measurement epoch. The following
section outlines an alternative approach to this problem using direct transcription methods.
3.2 Hermite-Simpson Polynomial Transcription
Explicit integration of the dynamics in (2) is often difficult to the unpredictable and chaotic behavior of cislunar
orbits over long time scales. Implicit integration methods have been shown to significantly reduce sensitivity to these
dynamics in IOD problems [6]. First, the domain is subdivided into N − 1 discrete intervals in time

tk = t1 < t2 < · · · < tN (18)

Similarly, the discretized states x k and controls u k are appended as unknown parameters. With Hermite-Simpson
collocation, the state is approximated by cubic polynomials [22], requiring the following constraints at each node

1  
gk xk , xk+1 , uk , uk+1/2 uk+1 = xk+1 − xk − hk f k + 4 f k+1/2 + f k+1 = 0 (19)
6
where f k = f (xxk , u k ,tk ) is the system dynamics evaluated at node k, and hk = tk+1 − tk . By construction, Hermite
polynomials require a value for the state and control at the midpoint of each interval. In compressed form, the value
of the state at the midpoint can be found by an interpolation of the trajectory at time tk+1/2 = 21 (tk + tk+1 ) as

1 hk
x k+1/2 = (xxk + x k+1 ) + ( f k − f k+1 ) (20)
2 8
However, since the control is interpolated with a quadratic spline, the controls at the midpoint of each interval must
included as decision variables.
The collocation equations in (19) enforce implicit integration at the boundaries and midpoint of each time interval.
However, the measurement residuals and admissible region constraints in (17) must be evaluated at each observation
epoch, which may not correspond to the discretization time of a node. This problem is accounted for in [6] by
constructing a spline approximation of the solution between nodes. First define τi = ti − tk as the time increment from
the nearest neighboring node. Using the Hermite-Simpson method, at measurement epoch an interpolating polynomial
can be constructed as
1   1  
x (τi ) = x k + f k τi + −3 f k + 4 f k+1/2 − f k+1 τi2 + 2 2 f k − 4 f k+1/2 + 2 f k+1 τi3 (21)
2hk 3hk

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
Similarly, a quadratic spline interpolation of the control variables gives

2 hk 4 2 hk
u (τi ) = u k (τi − )(τi − hk ) − 2 u k+1/2 τi (τi − hk ) + 2 u k+1 τi (τi − ) (22)
h2k 2 hk hk 2

These expressions allow for a solution to be evaluated between nodes. When evaluating the constraint Jacobian, the
derivatives must be included with respect to neighboring nodes k and k1 . See [6] for further details.
Before constructing the full NLP problem, we must also discretize integral quantities. Betts [23] describes two ap-
proaches to this problem. The first is to introduce an auxiliary state variable with rates of change (derivative) equal
to the functional of the integral in (16). The state is then integrated implicitly by the collocation equations. However,
this method increases the dimensionality unnecessarily. A more efficient approach is to utilize quadrature equations
to redefine the integral in terms of known quantities. For the Hermite-Simpson transcription, a quadrature rule can be
defined as
Z tf N−1
hk
u (τ)T u (τ) dτ ≈ ∑ (wk + 4wk+1/2 + wk+1 ) = b T q (23)
t0 k=1 6

For convenience we introduce the variable w(τ) = kuu(τ)k22 . In their factored forms, the coefficient vectors are defined
as
1
b T = h1 4h1+1/2 (h1 + h2 ) 4h2+1/2 (h2 + h3 ) . . . (hN−2 + hN−1 ) 4hN−1+1/2 hN−1

(24)
6
and  
w1
 w1+1/2 
 
 w2 
 
 w2+1/2 
 
q =  w3  (25)
 
 .. 

 . 

 wN−1 
 
wN−1+1/2 
wN
With these definitions in place, we have the necessary pieces to formulate a discrete approximation of the continuous
time optimization problem in (17). This problem can be approached using efficient interior point solvers such as
IPOPT [24] and SNOPT [25]. A major benefit of these methods is their large region of convergence and significantly
reduced sensitivity to initial guess errors when compared to direct shooting methods. Integration constraints are
satisfied using higher-order Hermite polynomials. The following section describes the NLP problem structure and
conditions for optimality.
3.3 Nonlinear Programming
The discretization of the problem in (17) is achieved using collocation equations (19) and interpolation (21)–(22) to
evaluate the solution between nodes. The constraints represent implicit integration of the system dynamics, as well as
measurement constraints for measurement association across UCTs. Finally, the integral control cost is represented
by the quadrature rule in (23).
The NLP problem structure is defined as follows

minimize b T q (uu1 , . . . , u N )




subject to: g k x k , x k+1 , u k , u k+1/2 , u k+1 = 0

N,m
Problem P (26)

 c L ≤ c (xxi ,ti ) ≤ cU
for i = 1, . . . , p, k = 1, . . . , N − 1


where the PN,m notation refers to an N-discretization with m-dimensional algebraic states. In order to better understand
this problem, it is helpful to develop the structure of the underlying sparsity. This information is commonly required
by NLP solvers to improve search direction computation. The NLP objective is only a function of the control variables,

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
Fig. 3: Illustration of the impulsive ∆V scenario for optical angles IOD. Spacecraft maneuvers are modeled by discrete
velocity increments at the collocation points.

giving the following derivatives

 
∂J
= 0T (27a)
∂ xk 1×n
 
∂J
= 2bk u Tk (27b)
∂ uk 1×m
 
∂J
= 2bk+1/2 u Tk+1/2 (27c)
∂ u k+1/2 1×m

where bk represents the k-th entry of the coefficient vector b in (24). It is convenient to define an augmented parameter
vector consisting of the total unknowns

X = (xx1 , x 2 , . . . , x N , u 1 , u 1+1/2 , u 2 , . . . , u N )T (28)

Let us also define an augmented constraint vector of the collocation and measurement constraints

X ) = (gg1 , g 2 , . . . , g N−1 , c 1 , c 2 , . . . , c p )T
z (X (29)

The constraint Jacobian has the following sparsity structure, where X’s represent blocks of non-zero entries.
 
XX XXX
 X X XXX 
 .. .. .. .. .. 
 . . . . . 
∂z  X X X X X
∼ X X ··· (30)
 

∂X  . .. 
 .
 . .


 ··· X X 
| {z }| {z }
x 1 ,xx2 ,...,xxN u 1 ,uuk+1/2 ,uu2 ,...,uuN

These expressions serve to illustrate the NLP gradient and Jacobian structures, which can be exploited to speed up
matrix inverse operations of the interior point solver. Problem (26) represents the NLP formulation of the IOD and
measurement association problem with maneuvering objects. The approach is developed for general nonlinear dynam-
ical systems, although results are focused on cislunar tracking problems in the following sections.
3.4 Impulsive Delta-V Formulation
The integrated control energy cost (16) can be used to solve IOD problems with low-thrust maneuvers. In practice,
spacecraft may operate with ∆V impulses that may not be well-suited to continuous thrust models. This problem can
be approached by reformulating the NLP structure for velocity inputs instead of a spline interpolation of the control,
which assumes the input is a smooth function. New variables are introduced to model impulsive increments in velocity
at each node. A high-level illustration of this problem is shown in Fig. 3.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 (a) Tanh smoothing (b) Slack variable formulation

Fig. 4: Comparison of absolute value smoothing methods and slack variable formulation. Each method avoids discon-
tinuities through the origin.

To begin, suppose the systems dynamics (11) can be stated as

ẋx(t) = f (xx(t),t) (31)

where control inputs u (t) are omitted. In order to ensure consistency with a transcription of the dynamics at each node,
the collocation constraints (19) are modified to include the velocity variable ∆vvk ∈ Rm as follows
1  T
g k (xxk , x k+1 , ∆vvk ) = x k+1 − x k − hk ( f k + 4 f k+1/2 + f k+1 ) + 0 I m ∆vvk = 0 (32)
6
In this case f k = f (xxk ,tk ) represents the dynamics evaluated without control input at node tk . Note that, in contrast to
the continuous control case, midpoint values of ∆vvk+1/2 are not required because the quadratic spline interpolation is
no longer used. This feature reduces the problem dimensionality by removing these decision variables from the NLP
formulation.
In order to properly model impulsive ∆V maneuvers, the control objective (16) must be reformulated to account for
total velocity increment, as the integral cost is based on an assumption for smooth functions. We instead consider a
sum of absolute values of velocity increments in each dimension, written as
N−1
J= ∑ k∆vvk k1 (33)
k=1

Note that this expression no longer relies on the quadrature relationships in the preceding section. One immediate
obstacle with this new objective lies in the implementation of the 1-norm over the control variables v k . Derivatives
are ill-defined at a value of zero, which is a common case over periods without maneuvers. One approach to this
problem is to approximate the absolute value operation with a smooth function, such as hyperbolic tangent or square-
root smoothing [22]. However, these methods are sensitive to smoothing parameters, making convergence difficult to
assess.
A more rigorous approach to this problem is to introduce slack variables that remove discontinuity from the objective.
This process also retains the structure of the original problem. The approach is outlined as follows. First, write the
(+) (−)
arguments of (33) in terms of new variables s k ≥ 0 and s k ≥ 0 as
m
(+) (−)
k∆vvk k1 = ∑ s k, j + s k, j (34)
j=1

We next enforce a constraint that relates the slack variables to the velocity variables
(+) (−)
s k − s k − ∆vvk = 0 (35)

The slack variables can be interpreted as modeling the left and right (negative and positive) sides of the absolute value
operation. By including these quantities as new decision variables, we have effectively removed discontinuities from
the objective and its gradient. An illustration of each approach is given in Fig. 4.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
The resulting NLP formulation with impulsive maneuvers is stated as follows
  
N−1 m (+) (−)


 minimize ∑ k=1 ∑ j=1 s k, j + s k, j

 subject to: g k (xxk , x k+1 , ∆vvk ) = 0

Problem PN,3m (+) (−)
s k − s k − ∆vvk = 0 (36)

c L ≤ c (xxi ,ti ) ≤ cU




for i = 1, . . . , p, k = 1, . . . , N − 1


The augmented decision variable for this problem includes the unknown states, velocity increments, and slack variables
(+) (+) (+) (−) (−) (−)
X = (xx1 , x 2 , . . . , x N , ∆vv1 , ∆vv2 . . . , ∆vvN−1 , s 1 , s 2 , . . . , s N−1 , s 1 , s 2 , . . . , s N−1 )T (37)

As in the preceding section, we can develop expressions for the objective gradients with respect to the unknowns at
each node
   
∂J T ∂J
=0 , = 0T (38a)
∂ x k 1×n ∂ ∆vvk 1×m
! !
∂J ∂J
(+)
= 1T, (−)
= 1T (38b)
∂ sk 1×m
∂ s k 1×m

These expressions clearly show the benefit of the slack variable approach in removing discontinuity of the 1-norm
objective. The sparsity structure of the Jacobian matrix can be developed analogous to (30), but is omitted from this
paper for brevity.
3.5 Measurement Uncertainty
Optical angles measurements are often subject to small errors from image processing and astrometry. It is worth noting
that the constraint in (14) does not account for these uncertainties, instead requiring that the measurement residual be
driven to exactly zero. This approach could be modified to account for measurement noise. Consider the measurement
residual ζ as a normal random variable with E[ζζ ] = 0 . Next, take the Cholesky decomposition of the inverse of the
measurement covariance to be R −1 = L T L . Then it follows that
T T
E[ζζ R−1 ζ ] = E[ζζ L T L ζ ]
Lζ k22
 
= E kL (39)

which implies a chi-square distribution with m degrees of freedom. This form can be related to a statistical significance
2 value and requiring that
level γ by finding the critical χγ,m
T
ζ R−1 ζ − χγ,m
2
≤0 (40)

This statistical relationship can be used to replace the equality constraint in (15). For brevity, however, the implemen-
tation of this approach is left for future work.

4. RESULTS

Examples of SSA-related applications of the collocation IOD algorithm developed in this paper are shown in the
following sections. These results generate a “truth” simulation of a maneuvering spacecraft in order to simulate
optical angles inputs for the IOD algorithm. No additional information about the target’s orbit or thrust capability is
assumed, including no initial guess for the state or control variables. The example scenarios are chosen to illustrate
the capability of the algorithm with continuous and impulsive thrust models.
4.1 Case 1: Low-Thrust Orbit Insertion to L4
The first example IOD scenario chosen for study begins with the target in an L4 long period orbit. The target actuates
a low-thrust maneuver to insertion at L4 over course of approximately 120 days. An observer placed in a periodic L1
Northern Halo orbit takes optical angles-only measurements of the target. Figure 5 illustrates the relative observing

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 Parameter Value
Initial Guess L4
Nodes 40
Free Variables 280
Constraints 255
Iterations 149
CPU Time 31 s
Objective 5.225E-04
NLP Error 9.684E-05

Table 1: NLP solution metrics with

Fig. 5: Low-thrust target orbit insertion to L4 (L1 Northern Halo observer). IPOPT solver.

geometry of this problem. Note that axis scaling is not equal, exaggerating the z dimension for clarity. The observer’s
orbit is plotted along with red markers indicating the location of measurement epochs along its orbit. The IOD output,
including reconstructed maneuvers, is described by its solution at the discretized collocation nodes. A spline interpo-
lation of this solution is shown in order to visualize a continuous solution of this trajectory. The true target trajectory
is plotted, although the estimates are close enough to make visualizing errors difficult. Table 1 gives parameters from
the solution of this IOD problem. The NLP solution indicates satisfactory optimality conditions, indicating the solver
has found a local minimum. A single-threaded solution takes approximately 30 seconds.
An observation history in terms of right ascension and declination are give in Fig. 6. Both the NLP solution nodes
and a spline interpolation of the solution are shown. It is emphasized that angles observations are taken only at
isolated measurement epochs (red markers). These plots illustrate the complex and varying dynamics of the relative
observing geometry over long time scales. For example, the oscillating trends in declination reflect the observer
transiting multiple periods of its orbit between observations of the target. In the process of finding a feasible solution,
the measurement constraints (15) must be precisely met. Therefore one may infer these observations are correlated by
the admissible region bounds.
An important output of the IOD algorithm is the reconstructed control history of a maneuvering target. Figure 7 shows
both true and estimated control magnitudes over time. The simulated truth trajectory actuates a total ∆V of 108 m/s.
In comparison, the IOD algorithm predicts a total ∆V of 98.5 m/s, corresponding to about a 9% error. This slight
under-estimation is readily explained by the control trends in Fig. 7. It is easy to see that both the magnitude and
shape of the control inputs are captured well by the algorithm. However, at the initial and final times, the control
estimate noticeably under-predicts the truth. This result is typical, as the true target trajectory is subject to initial and
final constraints (i.e., insertion to L4), whereas the IOD algorithm is free to optimize the endpoints. The reconstructed
trajectory does not precisely terminate at L4, instead finding a nearly ballistic trajectory producing the required angles
measurements.
This example scenario illustrates application of the collocation IOD algorithm with actively maneuvering targets.
Trajectory optimization with low-thrust models represents a wide range of targets. Extension of IOD methods for both
non-Keplerian dynamics and thrusting spacecraft provides significant SSA capability in cislunar space.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 Fig. 6: Optical angles measurements of reconstructed low-thrust trajectory.

Fig. 7: True and reconstructed low-thrust control inputs.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 Parameter Value
Initial Guess L2
Nodes 25
Free Variables 375
Constraints 243
Iterations 26
CPU Time 7.9 s
Objective 1.3780E-03
NLP Error 1.6722E-05

Table 2: NLP solution metrics with

Fig. 8: Impulsive ∆V station-keeping near L2 (L1 Northern Halo observer). IPOPT solver.

4.2 Case 2: L2 Southern Halo Station-Keeping

This example considers an impulsive ∆V model for spacecraft maneuvers within the collocation IOD algorithm. A
target spacecraft in an L2 Southern Halo orbit performs a series of ∆V maneuvers along its orbit. The maneuvers are
designed to correct an initial state deviation of
 T
δ x 0 = 1 km 0 0 1 m/s −1 m/s 0 (41)

over the course of one orbit period. An observer in an L1 Northern Halo orbit takes angles and angle rates mea-
surements of the target. A total of four measurements are taken, with the goal of correlating these observations and
reconstructing the ∆V control history.
An illustration of this example scenario is outlined in Fig. 8. Along with the observer’s orbit, a reconstruction of
the target orbit from the optical measurements is shown. The results show close agreement with the simulated truth
trajectory. Table 2 gives problem parameters and performance metrics from the NLP solver. In comparison to Table 1,
the impulsive problem has a greater number of free variables, even though fewer collocation nodes were created. This
effect is the result of augmenting additional slack variables and constraints in order to enforce the 1-norm objective
in (36). However, the solution still converges in fewer iterations, taking less than 10 seconds. The objective and NLP
error are below tolerances, indicating that the solver has converged to a local minimum. The reconstructed angles and
angle rates measurements in Fig. 9 show close agreement between the predicted and true optical measurements across
the trajectory.
The reconstructed ∆V solution is shown in Fig. 10. Results are presented as stem plots to emphasize the discrete
nature of quantities (values are no longer sampling a continuous curve). The truth simulation has two distinct peaks,
each corresponding to station-keeping maneuvers at the start of the orbit and a smaller correction about seven days
later. The reconstructed maneuver magnitudes are a close match to truth in terms of both magnitude and timing. The
predicted total ∆V is 1.371 m/s. When compared to the true total ∆V of 1.644 m/s, this corresponds to approximately
16% estimation error. The slack variable approach described in preceding sections replaces the 1-norm operation for
minimizing the sum of absolute ∆V ’s. The slack variables are zero everywhere except at indices where maneuvers are
detected. The sign of the maneuver is determined by either the positive or negative branch corresponding to each pair

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 (a) Angles (b) Angles rates

Fig. 9: Optical angles measurements of reconstructed station-keeping trajectory.

Fig. 10: True and reconstructed impulsive ∆V control inputs.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
 Fig. 11: Convergence qualities for the two example cases.

of slack variables. These results indicate the proposed method is successful in modeling impulsive maneuvers without
relying on approximate smoothing methods.
This example outlines application of the collocation IOD algorithm to a station-keeping scenario with impulsive con-
trol. Both the target orbit and its maneuver history are successfully reconstructed from angles and angle rates obser-
vations. In experimentation, it was found that angles-only observations were not sufficient for this problem without
increasing the number of measurements. This distinction is important to note when planning observations for poten-
tially maneuvering objects.

5. DISCUSSION AND FUTURE WORK

The example cases illustrate the capability of direct transcription methods for solving challenging IOD problems. The
formulation accounts for non-Keplerian motion, maneuver reconstruction, and observation association. While these
capabilities advance the state-of-the-art for cislunar SSA, they come with a moderate cost in computing overhead.
The results in [6] demonstrate convergence for non-maneuvering objects on the order of seconds. The inclusion of
additional control variables and constraints moderately increases the complexity of the problem and therefore solution
time. A comparison of the convergence rate for each example scenario in the preceding section is shown in Fig. 11. It
is interesting to note that, although the problem scale is larger in Case 2, the solution actually converges faster. This
is likely due to the effect of having a shorter timescale, as well as fewer total maneuvers compared to the continuous
thrust example in Case 1. The first example also appears to traverse into a local infeasible region, indicated by both
a “valley” effect coupled with a spike in the dual infeasibility. The solver is eventually able to escape this region and
converge on a local minimum. The second example does not exhibit the same behavior, leading to faster solution
times. These plots are useful to evaluate the overall progress of the solution, including both the constraint violation
and NLP error.
While the results of the collocation-based IOD approach are promising, the method relies on some moderate assump-
tions. First, the system dynamics are assumed known, and any unexplained motion is accounted for with control
inputs. In reality, dynamic mismodeling is likely to contribute to the system, such as four-body gravity effects or
solar radiation pressure. Using higher-fidelity models can alleviate this effect in part. Second, although the spacecraft
thrust capability is not assumed known, the distinction between impulsive and continuous thrust maneuvers must be
made prior to solution. Finally, measurement association may not be fully explained by admissible region constraints.
More robust methodologies using models of uncertainty are likely to provide higher confidence in UCT association
for cislunar objects.
The example cases serve to illustrate application of the collocation IOD algorithm to representative cislunar observing
scenarios. Further study is needed to determine the required measurement intervals for IOD and to understand the
effect of angles-only versus angles rates measurements. Future work should also consider observability metrics for
determining if follow-on observation of an IOD solution is required. Mass-optimal or fuel-limited trajectories could

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
also be studied. Finally, further investigation of the underlying NLP sparsity may provide a means for more efficient
numerical solution.

6. CONCLUSION

This paper presents a novel collocation-based IOD method for maneuvering objects in cislunar space. Problem def-
initions and models are provided. A discretization scheme using Hermite-Simpson polynomial transcription is then
described. The resulting NLP problem structure is discussed, and methods for solution are given. A modification
of the algorithm for determining impulsive velocity changes is also shown. Representative examples are solved for
cislunar observing applications with optical measurements, and results show close agreement with simulated quanti-
ties. The methodology developed in this paper enhances IOD capabilities in non-Keplerian systems. These findings
are expected to support spaceflight safety and operational awareness by advancing orbit determination capabilities in
cislunar space.

ACKNOWLEDGEMENTS

Support for this research was provided through an AFRL Cooperative Agreement (award No. FA9453-22-2-0050) at
the University of Colorado Boulder. The authors gratefully acknowledge technical discussion and advice from Dr. Jill
Bruer and Capt Lester Tuck at AFRL/RVSW.

REFERENCES

[1] Pedro Escobal. Methods of orbit determination. Methods of orbit determination, 1970.
[2] Samuel Wishnek, Marcus J Holzinger, and Patrick Handley. Robust cislunar initial orbit determination. In AMOS
Conf. Proc, 2021.
[3] Mark Bolden, I Hussein, H Borowski, R See, and E Griggs. Probabilistic initial orbit determination and object
tracking in cislunar space using optical sensors. In Advanced Maui Optical and Space Surveillance Technologies
(AMOS) Conference, pages 27–30, 2022.
[4] C Channing Chow, Charles J Wetterer, Jason Baldwin, Micah Dilley, Keric Hill, Paul Billings, and James Frith.
Cislunar orbit determination behavior: processing observations of periodic orbits with gaussian mixture model
estimation filters. The Journal of the Astronautical Sciences, 69(5):1477–1492, 2022.
[5] Sajjad Kazemi, Nasser L Azad, K Andrea Scott, Haroon B Oqab, and George B Dietrich. Orbit determination
for space situational awareness: A survey. Acta Astronautica, 2024.
[6] Casey R Heidrich and Marcus J Holzinger. Universal angles-only cislunar orbit determination using sparse collo-
cation. In Proceedings of the Advanced Maui Optical and Space Surveillance (AMOS) Technologies Conference,
page 17, 2023.
[7] Brynjulf Owren and Marino Zennaro. Derivation of efficient, continuous, explicit runge–kutta methods. SIAM
journal on scientific and statistical computing, 13(6):1488–1501, 1992.
[8] Keric Hill, Kyle Alfriend, and Chris Sabol. Covariance-based uncorrelated track association. In AIAA/AAS
Astrodynamics Specialist Conference and Exhibit, page 7211, 2008.
[9] Marcus J Holzinger, K Kim Luu, Chris Sabol, and Keric Hill. Uncorrelated-track classification, characterization,
and prioritization using admissible regions and bayesian inference. Journal of Guidance, Control, and Dynamics,
39(11):2469–2484, 2016.
[10] Johnny L Worthy III, Marcus J Holzinger, and Daniel J Scheeres. An optimization approach for observation asso-
ciation with systemic uncertainty applied to electro-optical systems. Advances in Space Research, 61(11):2709–
2724, 2018.
[11] Marcus J Holzinger, Daniel J Scheeres, and Kyle T Alfriend. Object correlation, maneuver detection, and char-
acterization using control distance metrics. Journal of Guidance, Control, and Dynamics, 35(4):1312–1325,
2012.
[12] Daniel P Lubey and Daniel J Scheeres. An optimal control based estimator for maneuver and natural dynam-
ics reconstruction. In Proceedings of the 2013 Advanced Maui Optical and Space Surveillance Technologies
Conference, 2013.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com
[13] Eusebius J Doedel, Volodymyr A Romanov, Randy C Paffenroth, Herbert B Keller, Donald J Dichmann, Jorge
Galán-Vioque, and André Vanderbauwhede. Elemental periodic orbits associated with the libration points in the
circular restricted 3-body problem. International Journal of Bifurcation and Chaos, 17(08):2625–2677, 2007.
[14] Johnny L Worthy III and Marcus J Holzinger. Incorporating uncertainty in admissible regions for uncorrelated
detections. Journal of Guidance, Control, and Dynamics, 38(9):1673–1689, 2015.
[15] Woosang Park and Kyle T Alfriend. Dynamic model fidelity effects on covariance based track association.
Journal of Guidance, Control, and Dynamics, pages 1–11, 2024.
[16] Kyle J DeMars and Moriba K Jah. Probabilistic initial orbit determination using gaussian mixture models.
Journal of Guidance, Control, and Dynamics, 36(5):1324–1335, 2013.
[17] Kohei Fujimoto and Daniel J Scheeres. Applications of the admissible region to space-based observations.
Advances in Space Research, 52(4):696–704, 2013.
[18] Mohammed A Ghazy and Brett Newman. Analytic theory for high-inclination orbits in the restricted three-body
problem. Journal of guidance, control, and dynamics, 33(2):565–583, 2010.
[19] Kenta Oshima, Stefano Campagnola, and Tomohiro Yanao. Global search for low-thrust transfers to the moon in
the planar circular restricted three-body problem. Celestial Mechanics and Dynamical Astronomy, 128:303–322,
2017.
[20] Bob Schutz, Byron Tapley, and George H Born. Statistical orbit determination. Elsevier, 2004.
[21] Andris D Jaunzemis, Midhun V Mathew, and Marcus J Holzinger. Control cost and mahalanobis distance binary
hypothesis testing for spacecraft maneuver detection. Journal of Guidance, Control, and Dynamics, 39(9):2058–
2072, 2016.
[22] Matthew Kelly. An introduction to trajectory optimization: How to do your own direct collocation. SIAM Review,
59(4):849–904, 2017.
[23] John T Betts. Practical methods for optimal control and estimation using nonlinear programming. SIAM, 2010.
[24] Lorenz T Biegler and Victor M Zavala. Large-scale nonlinear programming using ipopt: An integrating frame-
work for enterprise-wide dynamic optimization. Computers & Chemical Engineering, 33(3):575–582, 2009.
[25] Philip E Gill, Walter Murray, and Michael A Saunders. Snopt: An sqp algorithm for large-scale constrained
optimization. SIAM review, 47(1):99–131, 2005.

Copyright © 2024 Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS) – www.amostech.com