23-P12

# Autonomous Local Catalog Maintenance of Close Proximity Satellite Systems on
Closed Natural Motion Trajectories††thanks: Approved for public release;
distribution is unlimited. Public Affairs approval AFRL20230485

Christopher W. Hays Ph.D. Candidate, Aerospace Engineering Department, Embry-
Riddle Aeronautical University, 1 Aerospace Blvd., Daytona Beach, FL, 32114.
Kristina Miller Ph.D. Candidate, Electrical and Computer Engineering
Department, University of Illinois Urbana-Champaign, 1308 W Main St, Urbana,
IL 61801 Alexander Soderlund Research Aerospace Engineer, Space Vehicles
Directorate, Air Force Research Laboratory, 3550 Aberdeen Dr. SE, Kirtland
AFB, NM Sean Phillips Troy Henderson Research Mechanical Engineer, Space
Vehicles Directorate, Air Force Research Laboratory, 3550 Aberdeen Dr. SE,
Kirtland AFB, NMAssociate Professor, Aerospace Engineering Department, Embry-
Riddle Aeronautical University, 1 Aerospace Blvd., Daytona Beach, FL., 32114.

###### Abstract

To enable space mission sets like on-orbit servicing and manufacturing, agents
in close proximity maybe operating too close to yield resolved localization
solutions to operators from ground sensors. This leads to a requirement on the
systems need to maintain a catalog of their local neighborhood, however, this
may impose a large burden on each agent by requiring updating and maintenance
of this catalog at each node. To alleviate this burden, this paper considers
the case of a single satellite agent (a chief) updating a single catalog. More
specifically, we consider the case of numerous satellite deputy agents in a
local neighborhood of a chief, the goal of the chief satellite is to maintain
and update a catalog of all agents within this neighborhood through onboard
measurements. We consider the agents having relative translational and
attitude motion dynamics between the chief and deputy, with the chief centered
at the origin of the frame. We provide an end-to-end solution of the this
problem through providing both a supervisory control method coupled with a
Bayesian Filter that propagates the belief state and provides the catalog
solutions to the supervisor. The goal of the supervisory controller is to
determine which agent to look at and at which times while adhering to
constraints of the chief satellite. We provide a numerical validation to this
problem with three agents.

##  1 Introduction

Space is a key component to the global economy by facilitating trillions of
dollars, annually, in the global market due to GPS alone [1]. Due to the lower
barrier to entry, space is becoming increasingly congested and the current
space surveillance network may not be able to fully maintain the entire
catalog of small objects in the future [2]. Therefore, future satellites may
have to maintain a local catalog of nearby space objects within close
proximity to assure safety of the systems and autonomously decide if any
evasive actions would be required. Recently, readiness levels for space
trusted autonomy are studied.[3] In this paper, we consider the problem of
local area catalog maintenance for satellite systems. This work may be
particularly useful for such mission types as in space servicing and
manufacturing, wherein a “foreman satellite” is tasked with maintaining local
awareness of a formation of satellite robots.

In the literature, the topic of satellite localization has been thoroughly
studied for close proximity (relative motion) satellite systems [4, 5, 6, 7,
8, 9, 10], and ground-based sensor networks [11, 12]. Due to the number of
space objects compared to the number of ground based sensors, the tasking of
ground sensors must be intelligently chosen [13, 14]. Situational awareness,
in general, can be defined as “the perception of the elements in the
environment within a volume of time and space, the (organizational)
comprehension of their meaning, and the projection of their status in the near
future” [15]. In particular, space situational awareness is a growing topic of
interest [16]. Space situational awareness (a.k.a. space domain awareness, or
SDA) has been expanded to encompass all space environmental impacts as well,
which leads to a multidisciplinary domain of research that incorporates facets
of information fusion, collection-tasking and exploitation to better assess
anomaly identification and prediction.

In this article, we consider the case of a single controllable chief satellite
equipped with an electro-optical sensor and on-board estimating capabilities
that must track and maintain the states of several uncontrollable deputy
satellites, each of which lay on a closed, elliptical natural motion
trajectory [17] about the chief. The local catalog we wish to maintain is the
concatenated list of all the deputies’ positional and velocity states relative
to the chief. However, a deputy’s state is inherently uncertain. Generally,
these state uncertainties are represented through stochastic parameters, such
as a mean and state covariance matrix pair for a Gaussian probability
distribution, also known as the “belief state.” Through the use of a Bayesian
filter [18], the fusion between a deputy’s prior belief state and an
observation of that deputy results in a posterior belief state with reduced
uncertainty. Crucially, the chief cannot observe the states of all deputies
simultaneously and is only able to make measurements of deputies when they lie
within the chief’s field-of-view. Thus, the uncertainties of unobserved agents
will grow accordingly with the amount of time elapsed without measurement.
This engenders the core of the local catalog maintenance problem: how does the
chief know when and where to look in order to provide a catalog of deputy
states that are closest to their true values?

Further, this article provides a decision-making algorithm which supplies a
control trajectory to the chief satellite such that the uncertainties of all
tracked deputy agents are guaranteed to lie underneath a predefined bound. The
deputies are modeled to maneuver about the chief according to the well-known
Clohessey-Wiltshire dynamics, while the chief’s attitudinal dynamics are
modeled according to gas jet thruster equations. The state estimation is
performed using the extended Kalman filter, and the sequence of control inputs
to the chief are computed via Model Predictive Control (MPC). Finally, a
supervising decision logic framework is constructed that determines which
sequence of local agents to observe based on the estimated relative positions
and uncertainty. This decision logic framework is general enough to be paired
alongside a graph-based approach to a multi-agent version of a similar problem
presented by Hays et. al. [19].

The remainder of the paper continues with a section providing essential
mathematical notation and definitions, an overview of the model predictive
control process, and a description of the local catalog maintenance problem.
The succeeding section details the proposed autonomous supervisor solution to
the maintenance problem. The effectiveness of this method is then demonstrated
through numerical simulations of the local maintenance scenario, followed by
concluding remarks.

##  2 Preliminaries

###  2.1 Notations and Definitions

Let ${\mathbb{R}}$ denote the set of all real numbers and let $\mathbb{N}$
denote the set of natural numbers. A frame for an $n$-dimensional space is
defined by $n$ orthonormal vectors $\hat{\mathcal{I}}_{i}$, $i\in[1,n]$. In
this paper, we work with spaces where $n=3$. The Lie Group $\textrm{SO}(3)$ is
the set of all real invertible $3\times 3$ matrices that are orthogonal with
determinant $1$, in this work referred to as rotation matrices
$\mathcal{R}_{\to}\in\textrm{SO}(3):$

$\textrm{SO}(3):=\\{\mathcal{R}_{\to}\in{\mathbb{R}}^{3\times 3}\ |\
(\mathcal{R}_{\to})^{T}\mathcal{R}_{\to}=\mathbf{I}_{3\times 3},\
\det[\mathcal{R}_{\to}]=1\\}$

Given two coordinate frames $\mathcal{A}$ and $\mathcal{B}$, the rotation of
vector ${\mathbf{x}}^{\mathcal{A}}$ as measured in frame $\mathcal{A}$ to
${\mathbf{x}}^{\mathcal{B}}$ as measured in frame $\mathcal{B}$ is denoted by
$\mathcal{R}_{\mathcal{A}\to\mathcal{B}}$, that is,
${\mathbf{x}}^{\mathcal{B}}=\mathcal{R}_{\mathcal{A}\to\mathcal{B}}{\mathbf{x}}^{\mathcal{A}}$.
Given three frames $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$, the
angular velocity of $\mathcal{B}$ relative to $\mathcal{A}$ as measured in
frame $\mathcal{C}$ is denoted by
$\bm{\omega}_{\mathcal{A}\mathcal{B}}^{\mathcal{C}}$. For a vector
$\textbf{a}:=(a_{1}\ a_{2}\ a_{3})^{\top}\in{\mathbb{R}}^{3}$, the cross
product operator $\mathcal{S}({\cdot}):{\mathbb{R}}^{3}\to\mathfrak{so}(3)$:

$\mathcal{S}({\textbf{a}})=\begin{pmatrix}0&-a_{3}&a_{2}\\\ a_{3}&0&-a_{1}\\\
-a_{2}&a_{1}&0\end{pmatrix}$ (1)

If $\textbf{b}\in{\mathbb{R}}^{3}$, then the cross-product
$\textbf{c}=\textbf{a}\times\textbf{b}$ can be written
$\textbf{c}=\mathcal{S}({\textbf{a}})\textbf{b}=-\mathcal{S}({\textbf{b}})\textbf{a}$.
Here, the Lie algebra $\mathfrak{so}(3)$ is the set of all real $3\times 3$
skew-symmetric matrices and is the tangent space of $\textrm{SO}(3)$ at the
identity $\mathbf{I}_{3\times 3}$:

$\mathfrak{so}(3):=\\{\mathcal{S}({\textbf{a}})^{T}=-\mathcal{S}({\textbf{a}})\
|\ \textbf{a}\in{\mathbb{R}}^{3}\\}$

For a vector $\textbf{a}\in{\mathbb{R}}^{3}$, the $\mathcal{L}^{2}$ norm
operator is defined as $|\textbf{a}|:=\sqrt{\textbf{a}^{\top}\textbf{a}}$. A
multivariate Gaussian distribution of an $n$-dimensional random variable is
parameterized by a mean vector $\hat{{\mathbf{x}}}\in{\mathbb{R}}^{n}$ and a
symmetric, positive semidefinite covariance matrix
$\mathbf{P}\in{\mathbb{R}}^{n\times n}$. Given matrices
$A\in{\mathbb{R}}^{m\times n}$ and $B\in{\mathbb{R}}^{p\times q}$, the
Kroenecker product is denoted by $\otimes$ and is given by

$A\otimes B=\begin{pmatrix}a_{11}B&a_{12}B&\dots&a_{1n}\\\
a_{21}B&a_{22}B&\dots&a_{2n}\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{m1}B&a_{m2}B&\dots&a_{mn}B\end{pmatrix}.$

Given a vector $a\in{\mathbb{R}}^{m}$, the projection of $a$ onto
${\mathbb{R}}^{n}$, $n\leq m$ is denoted by $a\downarrow{\mathbb{R}}^{n}$.
Next, we specify a mapping between the vector ${\mathbf{x}}$ in Cartesian
space to the azimuth-elevation or (Az,El) space
$g({\mathbf{x}}):{\mathbb{R}}^{3}\to{\mathbb{S}}^{2}$

$g({\mathbf{x}})=\begin{bmatrix}\text{Az}\\\
\text{El}\end{bmatrix}=\begin{bmatrix}\tan^{-1}(\frac{s_{2}}{s_{1}})\\\
\sin^{-1}(s_{3})\end{bmatrix}$ (2)

where $s=\frac{{\mathbf{x}}}{|{\mathbf{x}}|}$ and ${\mathbb{S}}^{2}$ is the
unit sphere with angular coordinates $\text{Az}\in[-\pi,\pi]$ and
$\text{El}\in[-\frac{\pi}{2},\frac{\pi}{2})$ by astronomical convention.

###  2.2 Model Predictive Control

Consider the following control system with state $x_{k}\in{\mathbb{R}}^{n}$ at
timestep $k$ that evolves according to a discrete-time system with process
function $f:{\mathbb{R}}^{n}\times{\mathbb{R}}^{m}\to{\mathbb{R}}^{n}$ and
measurement function
$c:{\mathbb{R}}^{n}\times{\mathbb{R}}^{m}\to{\mathbb{R}}^{\ell}$ over
timesteps $\\{k,k+1,k+2,\dots\\}$, $k\geq 0$ as

$\displaystyle x_{k+1}$ $\displaystyle=f(x_{k},u_{k})$ (3) $\displaystyle
y_{k}$ $\displaystyle=c(x_{k},u_{k})$ (4)

where $u_{k}\in{\mathbb{R}}^{m}$ is the input applied to the plant while
$y_{k}\in{\mathbb{R}}^{\ell}$, $\ell\leq n$, is the output of the system.
Without loss of generality, we assume $f(0,0)=0$ and $c(0,0)=0$ for all $k\geq
0$. Let $N\in\mathbb{N}$, $N\geq 2$ be the finite horizon length. We define a
control policy $\textbf{u}_{k}:=\\{u_{0},u_{1},\dots,u_{i},\dots,u_{N-1}\\}$
at time $k$ with an additional index $i=0,\dots,N-1$. The output trajectory
resulting from $\textbf{u}_{k}$ according to (3)-(4) is
$\textbf{y}_{k}:=\\{y_{0},y_{1},\dots,y_{i},\dots,y_{N-1}\\}$. Let
$\textbf{y}_{k}(i)$ and $\textbf{u}_{k}(i)$ denote particular elements of
$\textbf{y}_{k}$ and $\textbf{u}_{k}$ at index $i$ given an initial condition
$x_{k}$, respectively. Let there be some nonnegative function
$h:{\mathbb{R}}^{\ell}\times{\mathbb{R}}^{m}\to{\mathbb{R}}_{\geq 0}$, known
as a stage cost, and a control objective function
$J:{\mathbb{R}}^{\ell}\times{\mathbb{R}}^{m}\to{\mathbb{R}}_{\geq 0}$ defined
as

$J(\textbf{y}_{k},\textbf{u}_{k}):=\sum^{N-1}_{i=0}h(\textbf{y}_{k}(i),\textbf{u}_{k}(i)).$
(5)

The goal of a model predictive controller is to, at timestep $k$, generate an
optimal control sequence
$\textbf{u}^{*}_{k}:=\\{u^{*}_{0},u^{*}_{1},\dots,u^{*}_{i},\dots,u^{*}_{N-1}\\}$
such that the objective function $J$ from (5) is minimized. The first input in
this sequence, $\textbf{u}^{*}_{k}(0)$, is applied to the system. At timestep
$k+1$, this minimization process is repeated and $\textbf{u}^{*}_{k+1}(0)$ is
applied. This is formally written as

$\displaystyle\begin{split}\underset{\textbf{u}_{k}}{\textnormal{minimize}}&\
J(\textbf{y}_{k},\textbf{u}_{k})\\\ \textnormal{subject to}&\ x_{0}=x_{k}\\\
&\ x_{i+1}=f(x_{i},u_{i}),\quad i=0,\dots,N-1\\\ &\ y_{i}=c(x_{i},u_{i}),\quad
i=0,\dots,N-1\\\ &\ (x_{i},u_{i})\in\mathscr{X}\times\mathscr{U},\quad
i=0,\dots,N-1\\\ \end{split}$ (6)

where $\mathscr{X}\subseteq{\mathbb{R}}^{n}$ is the closed set of state
constraints and $\mathscr{U}$ is the closed set of input constraints.
Generally, the stage cost function $h$ is chosen to represent a distance of
the output from some predefined desired output trajectory, i.e.
$\textbf{y}^{*}_{k}$.

##  3 Problem Statement

###  3.1 The Local Catalog Maintenance Problem

Consider a chief spacecraft with state space
$\mathcal{X}_{0}\subseteq{\mathbb{R}}^{n}$ and input space
$\mathcal{U}\subseteq{\mathbb{R}}^{m}$. The state
$\mathbf{x}_{0}\in\mathcal{X}_{0}$ evolves according to
$\dot{\mathbf{x}}_{0}=f(\mathbf{x}_{0},u)$ where $u\in\mathcal{U}$ is the
chief’s input, and $f$ is the chief’s dynamic function. Additionally, consider
some constraint function
$h:\mathcal{X}_{0}\times\mathcal{U}\to{\mathbb{R}}^{r}$. A state
$\mathbf{x}_{0}\in\mathcal{X}_{0}$ and input $u\in\mathcal{U}$ satisfies the
constraints when $h_{i}(\mathbf{x}_{0},\mathcal{U})\leq 0$ for each row
$i\in\\{1,\dots,r\\}$ in the constraint vector. Given an initial state
$\theta_{0}\in\mathcal{X}_{0}$ and an input function $u:\mathbb{R}_{\geq
0}\to\mathcal{U}$, a trajectory of the chief is given by
$\xi_{u,\theta_{0}}:\mathbb{R}_{\geq 0}\to\mathcal{X}$ and satisfies
$\xi_{u,\theta_{0}}(0)=\theta_{0}$ and
$\frac{d}{dt}\xi_{u,\theta_{0}}(t)=f(\xi_{u,\theta_{0}}(t),u(t))$.

The chief is tasked with maintaining a catalog of $d$ deputies. Each deputy
$i\in\\{1,\dots,d\\}$ has a state space
$\mathcal{X}_{i}\subseteq{\mathbb{R}}^{p}$. A deputy’s state
${\mathbf{x}}_{i}\in\mathcal{X}_{i}$ evolves according to
$\dot{{\mathbf{x}}}_{i}={\mathbf{A}}{\mathbf{x}}_{i}$ where
${\mathbf{A}}\in{\mathbb{R}}^{p\times p}$ is the state matrix. The state of
all deputies, or system state, is given by
${\mathbf{x}}:=[{\mathbf{x}}_{i}^{\top}\ \dots\
{\mathbf{x}}_{d}^{\top}]^{\top}$ and evolves according to
${\mathbf{x}}=(\mathbf{I}_{d\times d}\otimes{\mathbf{A}}){\mathbf{x}}$.

The chief is equipped with a recursive state estimator to track $d$ deputy
belief states, parameterized by concatenated mean
$\hat{{\mathbf{x}}}:=[{\mathbf{x}}^{\top}_{1}\ \cdots\
{\mathbf{x}}^{\top}_{d}]^{\top}$ and covariance block matrix
$\mathbf{P}:=\text{diag}\\{\mathbf{P}_{1},\ldots,\mathbf{P}_{d}\\}$. These
constructs evolve via their assumed dynamics models according to

$\displaystyle\dot{\hat{{\mathbf{x}}}}=(\mathbf{I}_{d\times
d}\otimes{\mathbf{A}})\hat{{\mathbf{x}}}$
$\displaystyle\dot{\mathbf{P}}=(\mathbf{I}_{d\times
d}\otimes{\mathbf{A}})\mathbf{P}(\mathbf{I}_{d\times
d}\otimes{\mathbf{A}})^{\top}+\mathbf{G}\mathbf{Q}\mathbf{G}^{\top},$ (7)

where, $\mathbf{Q}\in{\mathbb{R}}^{pd\times pd}$ is the process noise matrix
and $\mathbf{G}\in{\mathbb{R}}^{pd\times pd}$ is the noise input matrix. We
assume $\mathbf{G}=\mathbf{I}_{pd\times pd}$ for simplicity.

The chief is equipped with a sensor parameterized by its pointing vector in
the chief’s body frame $\mathbf{p}^{\mathcal{B}}$ and its angle of view
$\alpha$. The sensor can observe the states of the deputies within some
sensing field-of-view (FOV) dependent on $\mathbf{p}^{\mathcal{B}}$, $\alpha$,
and the chief’s state $\mathbf{x}_{0}$. The linear measurement model of the
sensor is given by $\mathbf{y}=C({\mathbf{x}},\mathbf{x}_{0}){\mathbf{x}}$
where $\mathbf{y}$ is the observed state of deputies, and $C$ is the
observation matrix dependent on the system state ${\mathbf{x}}$ and the
chief’s state $\mathbf{x}_{0}$. The observation matrix $C$ is a block matrix
given as

$C({\mathbf{x}},\mathbf{x}_{0})=\begin{bmatrix}C_{1}({\mathbf{x}}_{1},\mathbf{x}_{0})&\mathbf{\emptyset}_{p\times
p}&\cdots&\mathbf{\emptyset}_{p\times p}\\\ \mathbf{\emptyset}_{p\times
p}&C_{2}({\mathbf{x}}_{2},\mathbf{x}_{0})&\cdots&\mathbf{\emptyset}_{p\times
p}\\\ \vdots&\vdots&\ddots&\vdots\\\ \mathbf{\emptyset}_{p\times
p}&\mathbf{\emptyset}_{p\times
p}&\cdots&C_{d}({\mathbf{x}}_{d},\mathbf{x}_{0})\end{bmatrix},$ (8)

where each $C_{i}({\mathbf{x}}_{i},\mathbf{x}_{0})$ is given as

$C_{i}({\mathbf{x}}_{i},\mathbf{x}_{0})=\begin{cases}\mathbf{I}_{p\times
p}&\textrm{if ${\mathbf{x}}_{i}$ is observed}\\\ \mathbf{\emptyset}_{p\times
p}&\textrm{else.}\end{cases}$ (9)

It should be noted here that $\mathbf{\emptyset}_{p\times p}$ is an overloaded
operator. In the case that
$C_{i}({\mathbf{x}}_{i},\mathcal{R}_{\mathcal{B}\to\mathcal{I}})=\mathbf{I}_{p\times
p}$, $\mathbf{\emptyset}_{p\times p}$ takes on the value of
$\mathbf{0}_{p\times p}$, and in the case when
$C_{i}({\mathbf{x}}_{i},\mathcal{R}_{\mathcal{B}\to\mathcal{I}})=\mathbf{\emptyset}_{p\times
p}$, $\mathbf{\emptyset}_{p\times p}$ is the null set taking the form of an
empty matrix.

Whenever the chief makes an observation of a deputy $i\in\\{1,\dots,d\\}$ at a
time instant $k$, the measurement data is incorporated by the on-board
estimator to impulsively update the belief state such that

$(\hat{{\mathbf{x}}}_{i}(k^{-}),\mathbf{P}_{i}(k^{-}))\rightarrow(\hat{{\mathbf{x}}}_{i}(k^{+}),\mathbf{P}_{i}(k^{+})).$

As a result, the size of the covariance matrices of observed deputies will
necessarily decrease $f[\mathbf{P}_{i}(k^{+})]\leq f[\mathbf{P}_{i}(k^{-})]$
and the updated belief state will be recursively supplied to the propagator
until the next available measurement update. Here $f[\cdot]$ is any metric
quantifying the size of the covariance matrix or its representative
hyperellipsoid. However, because the sensor has a limited FOV, it may not be
possible for the chief to observe every deputy at once. Therefore, the chief
must maneuver to view each deputy to maintain a level of certainty in the
state estimate, that is $f[\mathbf{P}_{i}(k)]\leq\epsilon$ for all $k\geq 0$,
where $\epsilon$ is some positive value determined by the sensor FOV.

Then, the local catalog maintenance problem is the following. Consider a chief
satellite with dynamic function $f$ and constraint function $h$, a sensor
parametrized by pointing vector $\mathbf{p}^{\mathcal{B}}$ and area of view
$\alpha$. Given $d$ deputies with initial estimate
$\hat{\theta}\in\mathcal{X}$ and initial covariances
$\mathbf{P}^{0}\in{\mathbb{R}}^{pd\times pd}$, and initial chief state
$\theta_{0}\in\mathcal{X}_{0}$, find a control function
$c:\mathcal{X}_{0}\times\mathcal{X}\times{\mathbb{R}}^{pd\times
pd}\times\mathbb{R}_{\geq 0}\to\mathcal{U}$ such that for all $k\geq 0$,
$h_{i}(\xi_{c,\theta_{0}}(k),c(\xi_{c,\theta_{0}}(t)))\leq 0$ for all
$i\in\\{1,\dots,r\\}$, and $f[\mathbf{P}_{j}(k)]\leq\epsilon$.

#### Chief attitude dynamics

Consider a chief spacecraft in a circular orbit centered on the earth. While
it does not have thrusters for translational motion, it can change its
orientation using gas jet thrusters that produce external torques about the
chief’s center of mass. The chief’s attitudinal state is given by the Euler
angle set $\mathbf{x}_{0}:=[\psi,\ \theta,\ \phi]^{\top}$, respectively
denoted yaw, pitch, and roll, that parameterizes the rotation matrix from the
body frame to the inertial frame
$\mathcal{R}_{\mathcal{B}\to\mathcal{I}}=\mathcal{R}_{\mathcal{B}\to\mathcal{I}}^{\phi}\mathcal{R}_{\mathcal{B}\to\mathcal{I}}^{\theta}\mathcal{R}_{\mathcal{B}\to\mathcal{I}}^{\psi}$
using the 3-2-1 rotation sequence. The analogous angular velocities are
$\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}:=[\omega_{1},\
\omega_{2},\ \omega_{3}]^{\top}$. The chief inputs are given by the torques in
the body frame $u^{\mathcal{B}}:=[u_{1},\ u_{2},\ u_{3}]^{\top}$. The rotation
kinematics for the chief’s body-fixed frame $\mathcal{B}$ is given by

$\dot{\mathcal{R}}_{\mathcal{B}\to\mathcal{I}}=\mathcal{R}_{\mathcal{B}\to\mathcal{I}}\mathcal{S}({\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}}).$
(10)

The Euler angle set (locally) evolves according to

$\displaystyle\dot{\Gamma}$
$\displaystyle=\begin{pmatrix}-\cos(\psi)\tan(\theta)&-\sin(\psi)\tan(\theta)&-1\\\
\sin(\psi)&-\cos(\psi)&0\\\
-\cos(\psi)\sec(\theta)&-\sin(\psi)\sec(\theta)&0\\\
\end{pmatrix}\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}.$ (11)

Assuming that the chief spacecraft can be modeled as a rigid body and has an
inertia matrix ${\mathcal{J}}$ as measured in the body-fixed frame, the
application of Euler’s second law of motion gives the attitude dynamics

$\dot{\bm{\omega}}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}=-[{\mathcal{J}}]^{-1}\mathcal{S}({\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}})[{\mathcal{J}}\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}]+[{\mathcal{J}}]^{-1}u^{\mathcal{B}}.$
(12)

The angular velocity of the $\mathcal{B}$ frame relative to Hill’s frame
$\mathcal{H}$ using the property of relative angular velocity summation is

$\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{I}}=\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{I}}-\bm{\omega}_{\mathcal{I}\mathcal{H}}^{\mathcal{I}}$

with the time rate derivatives given as

$\dot{\bm{\omega}}_{\mathcal{H}\mathcal{B}}^{\mathcal{I}}=\dot{\bm{\omega}}_{\mathcal{I}\mathcal{B}}^{\mathcal{I}}-\dot{\bm{\omega}}_{\mathcal{I}\mathcal{H}}^{\mathcal{I}}.$
(13)

Using (10) and
$\bm{\omega}_{\mathcal{A}\mathcal{B}}^{\mathcal{B}}=\mathcal{R}_{\mathcal{A}\to\mathcal{B}}\bm{\omega}_{\mathcal{A}\mathcal{B}}^{\mathcal{A}}$
a manipulation of (13) yields

$\dot{\bm{\omega}}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}=\dot{\bm{\omega}}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}-\mathcal{R}_{\mathcal{H}\to\mathcal{B}}\dot{\bm{\omega}}_{\mathcal{I}\mathcal{H}}^{\mathcal{H}}-\mathcal{S}({\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}})\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}.$
(14)

Note that from the circular chief orbit assumption we have
$\dot{\bm{\omega}}_{\mathcal{I}\mathcal{H}}^{\mathcal{H}}=\mathbf{0}$ and that
$\bm{\omega}_{\mathcal{I}\mathcal{H}}^{\mathcal{H}}=(0\ 0\ \eta)^{\top}$ is
the constant angular velocity of the Hill’s frame relative to the inertial
frame where $\eta$ is the constant mean motion parameter. From the assumption
that the chief-fixed inertia matrix is aligned with its principal axes we have

${\mathcal{J}}=\begin{pmatrix}J_{1}&0&0\\\ 0&J_{2}&0\\\
0&0&J_{3}\end{pmatrix}.$

The combination of (12) and (14) yields the angular velocity dynamics as

$\dot{\bm{\omega}}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}=-[{\mathcal{J}}]^{-1}\mathcal{S}({\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}})[{\mathcal{J}}\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}]+[{\mathcal{J}}]^{-1}u^{\mathcal{B}}-\mathcal{S}({\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}})\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}.$

Using the relation
$\bm{\omega}_{\mathcal{I}\mathcal{B}}^{\mathcal{B}}=\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}+\bm{\omega}_{\mathcal{I}\mathcal{H}}^{\mathcal{B}}$
and cross product properties the above expression becomes

$\displaystyle\begin{split}\dot{\bm{\omega}}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}&=-[{\mathcal{J}}]^{-1}\big{(}\mathcal{S}({\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}})[{\mathcal{J}}\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}]+\mathcal{S}({\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}})[{\mathcal{J}}\bm{\omega}_{\mathcal{I}\mathcal{H}}^{\mathcal{B}}]+\\\
&\qquad\qquad\qquad\qquad\mathcal{S}({\bm{\omega}_{\mathcal{I}\mathcal{H}}^{\mathcal{B}}})[{\mathcal{J}}\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}]+\mathcal{S}({\bm{\omega}_{\mathcal{I}\mathcal{H}}^{\mathcal{B}}})[{\mathcal{J}}\bm{\omega}_{\mathcal{I}\mathcal{H}}^{\mathcal{B}}]\big{)}+\mathcal{S}({\bm{\omega}_{\mathcal{H}\mathcal{B}}^{\mathcal{B}}})\bm{\omega}_{\mathcal{I}\mathcal{H}}^{\mathcal{B}}+[{\mathcal{J}}]^{-1}u^{\mathcal{B}}.\end{split}$
(15)

#### Chief constraints

The chief also has some constraints it must follow while performing the
inspection. Typically, the torque inputs cannot exceed some $u_{\max}$, and
the angular velocity must be maintained under some $\omega_{\max}$. To
maintain safety, the sensor cannot point in the direction of the (assumed
static) sun vector $\mathbf{s}$ given in $\mathcal{I}$. Recall that
$\mathcal{R}_{\mathcal{B}\to\mathcal{I}}$ denotes the rotation matrix from the
chief’s body-fixed frame $\mathcal{B}$ to the inertial frame $\mathcal{I}$.
Then, the constraint function is given as

$h(\mathbf{x}_{0},u)=\begin{pmatrix}h_{1}(\mathbf{x}_{0},u)\\\
h_{2}(\mathbf{x}_{0},u)\\\ h_{3}(\mathbf{x}_{0},u)\\\
h_{4}(\mathbf{x}_{0},u)\\\ h_{5}(\mathbf{x}_{0},u)\\\
h_{6}(\mathbf{x}_{0},u)\\\ h_{7}(\mathbf{x}_{0},u)\\\
\end{pmatrix}=\begin{pmatrix}\mathbf{s}^{\top}(\mathcal{R}_{\mathcal{B}\to\mathcal{I}}\mathbf{p}_{\mathcal{B}})-\cos(\alpha)\\\
\|\omega_{1}\|-\omega_{\max}\\\ \|\omega_{2}\|-\omega_{\max}\\\
\|\omega_{3}\|-\omega_{\max}\\\ \|u_{1}\|-u_{\max}\\\ \|u_{2}\|-u_{\max}\\\
\|u_{2}\|-u_{\max}\\\ \end{pmatrix}.$ (16)

If $h_{1}(\mathbf{x}_{0},u)\leq 0$, then the angle between $\mathbf{s}$ and
$\mathbf{p}^{\mathcal{I}}=\mathcal{R}_{\mathcal{B}\to\mathcal{I}}\mathbf{p}^{\mathcal{B}}$
is at least $\alpha$, meaning the sun vector is not pointing in the sensor
FOV. If $h_{i}(\mathbf{x}_{0},u)\leq 0$, $i\in\\{2,3,4\\}$, then
$\|\omega_{j}\|\leq\omega_{\max}$ for $j\in\\{1,2,3\\}$. If
$h_{i}(\mathbf{x}_{0},u)\leq 0$, $i\in\\{5,6,7\\}$, then $\|u_{j}\|\leq
u_{\max}$ for $j\in\\{1,2,3\\}$.

#### Deputy relative motion dynamics

The state of the deputies are described in $\mathcal{H}$, and the
translational dynamics of the deputies are described by the Clohessy-Wiltshire
equations. The state of a deputy $i\in\\{1,\dots,d\\}$ is given by
${\mathbf{x}}_{i}=[r^{i}_{x}\ r^{i}_{y}\ r^{i}_{z}\ \dot{r}^{i}_{x}\
\dot{r}^{i}_{y}\ \dot{r}^{i}_{z}]^{\top}$, where $r^{i}_{x}$, $r^{i}_{y}$, and
$r^{i}_{z}$ are the $x$-, $y$-, $z$-positions of deputy $i$ in $\mathcal{H}$
respectively. The Clohessy-Wiltshire equations can be written in state-space
form as in (17).

$\dot{{\mathbf{x}}}_{i}=\begin{bmatrix}0&0&0&1&0&0\\\ 0&0&0&0&1&0\\\
0&0&0&0&0&1\\\ 3\eta^{2}&0&0&0&2\eta&0\\\ 0&0&0&2\eta&0&0\\\
0&0&-\eta^{2}&0&0&0\end{bmatrix}{\mathbf{x}}_{i}$ (17)

Each deputy’s relative state trajectory ${\mathbf{x}}_{i}$ is uncontrolled and
evolves under the following assumptions: (i) $\dot{r}^{i}_{y}=-2\eta
r^{i}_{x}$ and (ii) $\dot{r}^{i}_{x}=\frac{1}{2}\eta r^{i}_{y}$. Condition (i)
ensures that each deputy is in a periodic, closed natural motion trajectory
(NMT), and condition (ii) ensures that the NMT is a closed ellipse centered at
the origin of $\mathcal{H}$.

#### Observation condition

Recall from before that the deputy belief states are updated whenever an
observation is made. Given a deputy state ${\mathbf{x}}_{i}$, chief state
$\mathbf{x}_{0}$, with a sensor characterized by the pointing vector in the
body frame $\mathbf{p}^{\mathcal{B}}$ and the viewing angle $\alpha$. An
observation occurs if

$\cos^{-1}(\mathcal{R}_{\mathcal{B}\to\mathcal{H}}\mathbf{p}^{\mathcal{B}}\cdot({\mathbf{x}}_{i}\downarrow\mathcal{H}))\leq\alpha$

where
$\mathcal{R}_{\mathcal{B}\to\mathcal{H}}=\mathcal{R}_{\mathcal{I}\to\mathcal{H}}\mathcal{R}_{\mathcal{B}\to\mathcal{I}}^{\top}.$

##  4 Methodology

To solve the local catalog maintenance problem, this work proposes a multi-
layered closed-loop control strategy (see Fig.1) that incorporates continuous
dynamics with impulsive zero-order hold controls applied at discrete time
instances. That is, at a given time instant $t_{k}$, the state estimator fuses
predicted deputy belief states
$(\mathbf{\hat{x}}^{-}(t_{k}),\mathbf{P}^{-}(t_{k}))$ with available
observation data $y(t_{k})$ to yield posterior deputy belief states
$(\mathbf{\hat{x}}^{+}(t_{k}),\mathbf{P}^{+}(t_{k}))$ through a Bayesian
update scheme such as the commonly-used extended Kalman filter or unscented
Kalman filter.

Figure 1: Visual depiction of the proposed catalog maintenance operation. Note
that the sensor component is embedded within the plant component. The blue
dashed line indicates the recursive nature of the Bayesian filter, where
updated belief states are supplied to the propagation block for the next time
step.

###  4.1 Supervisor Algorithm

These updated beliefs are fed into a higher-level autonomous “supervisor”
which accounts for (i) deputy belief state information gaps (ii) prespecified
constraints and (iii) the current angular states of the chief. The supervisor
amalgamates this information and provides decision-making capabilities to a
lower-level controller, which drives the chief orientation trajectory by
administering torque inputs at discrete timesteps. Recall the local catalog
maintenance problem where a chief, with a sensor parametrized by
$\mathbf{p}^{\mathcal{B}}$ and $\alpha$, must maintain state estimates of $d$
deputies. The objective of the supervisor algorithm is to select the “optimal”
target deputy to track based on some cost function notated by $\mathcal{J}$. A
detailed breakdown of the supervisor algorithm can be seen in Algorithm 1.
Given the state estimates-covariance pairs of the deputies
$\\{(\hat{{\mathbf{x}}}_{i},\mathbf{P}_{i})\\}_{i=1}^{d}$ and the current
deputy $j$. Additionally, we introduce a new variable $\tau$ which “remembers”
the time that the supervisor switched to deputy $j$. The output is the deputy
$j^{\prime}$ to be viewed as well as $\tau^{\prime}$. Then, the supervisor
algorithm is the following. If some $\Delta$ time has not passed since the
last switching time, or the entropy of the current deputy has not fallen below
some threshold $\epsilon$, the current deputy $j$ and current $\tau$ is
returned. If not, then the entropy of all the deputies is computed. The deputy
to be viewed is chosen to be the one with maximum entropy, and $\tau$ is
updated to the current time. The supervisor cost function in this work is
considered to be the Shannon entropy of the covariance matrix, which is
assumed to be Gaussian; other metrics such as measures of uncertainty relating
to the deputies’ belief states, such as the $\det$ or $\mathrm{trace}$ of the
covariance matrix, or a measure of information gain post-fusion, such as the
posterior Fisher information [20] may also be utilized. The cost value for
each deputy is calculated, and then the deputy with the highest cost value is
selected to be the target deputy $j$. However, to prevent the possibility of
instantaneous switching, a form of hysteresis is applied. A target switch is
only triggered when the cost of the current target deputy falls below some
threshold, the potential other target deputy crosses over the same threshold,
and two target switches cannot occur within a time threshold,
$\delta\geq\gamma$, of each other, where $\delta$ is the time differential
between switches, and $\gamma$ is the desired time hysteresis threshold.

Input: $\\{(\hat{{\mathbf{x}}}_{i},\mathbf{P}_{i})\\}_{i=1}^{d}$, $j$, $\tau$

Output: $j^{\prime}$, $\tau^{\prime}$

1

2$\leftarrow\textsf{CurrentTime}$

3 if  _$(t-\tau\leq\Delta)\ \vee\
(\frac{k}{2}(1+\log(2\pi))+\log(|\mathbf{P}_{i}|) >\epsilon)$_ then

4 $j^{\prime}\leftarrow j$

5 $\tau^{\prime}\leftarrow\tau$

6

7else

8 for  _$i=1,\dots,d$_ do

9 $\mathcal{J}_{i}\leftarrow\frac{k}{2}(1+\log(2\pi))+\log(|\mathbf{P}_{i}|)$

10 $j^{\prime}\leftarrow\arg\max_{i\in\\{1,\dots,d\\}}\mathcal{J}_{i}$

11 $\tau^{\prime}\leftarrow t$

12

return $j^{\prime}$, $\tau^{\prime}$

Algorithm 1 Supervisor algorithm

###  4.2 Angle-Driven MPC

We implement a state tracking version of the model predictive controller
outlined in (6). Here we let $c(x_{k},u_{k}):=x_{k}$. Let there be two
positive definite square matrices $W_{1}\in{\mathbb{R}}^{2\times 2}$ and
$W_{2}\in{\mathbb{R}}^{3\times 3}$ and an Euler angle vector $z:=[\psi,\
\theta]^{\top}$. Note that the angle pair $\psi$ and $\theta$ were chosen due
to the mapping from the azimuth-elevation space. The sequence over discrete
steps from instance $k$ is denoted as
$\textbf{z}_{k}:=\\{z_{0},z_{1},\dots,z_{i},\dots,z_{N-1}\\}$ with a reference
angle trajectory (i.e. angle sequence associated with a current designated
target) is denoted as
$\textbf{z}^{r}_{k}:=\\{z^{r}_{0},z^{r}_{1},\dots,z^{r}_{i},\dots,z^{r}_{N-1}\\}$.
The MPC problem is defined with a quadratic stage cost as

$\displaystyle\begin{split}\underset{\textbf{u}_{k}}{\textnormal{minimize}}&\
J(\textbf{z}_{k},\textbf{u}_{k}):=\sum_{i=0}^{N-1}(\textbf{z}_{k}(i)-\textbf{z}^{r}_{k}(i))^{\top}W_{1}(\textbf{z}_{k}(i)-\textbf{z}^{r}_{k}(i))+(\textbf{u}_{k}(i))^{\top}W_{2}(\textbf{u}_{k}(i))\\\
\textnormal{subject to}&\ z_{0}=z_{k}\\\ &\ z_{i+1}=f_{z}(z_{i},u_{i}),\quad
i=0,\dots,N-1\\\ &\ (z_{i},u_{i})\in\mathscr{X}\times\mathscr{U},\quad
i=0,\dots,N-1.\\\ \end{split}$ (18)

Here, $f_{z}$ is equivalent to the dynamics expressed in (11) and (15).

##  5 Results

To demonstrate the efficacy of the proposed methodology, results of two
example scenarios are presented - a 3 deputy case and a 10 deputy case. For
this sample, the target deputy is selected via the Supervisor Algorithm
detailed in Algorithm 1. The MPC algorithm then uses the azimuth-elevation
tracks of deputy $j^{\prime}$ as the reference trajectory. The CasADi
optimization toolbox is used to solve the nonlinear optimal control problem
[21] using a direct single shooting method. The MPC cost function weighting
matrices were set as $\mathrm{W}_{1}=\mathbf{I}_{2\times 2}$ and
$\mathrm{W}_{2}=\mathbf{I}_{3\times 3}$. The torque constraint was set to be
$u_{max}=2\pi\ \text{Nm}$, the angular velocity constraint was set as
$\omega_{max}=\pi\ \text{rad/s}$ The example simulation consists of three
deputy agents following elliptical NMTs about the chief. The chief’s initial
attitude and angular velocity were randomly generated with the initial angular
velocity being restricted to a bounded box of $\pi/2\ \text{m/s}$ about the
origin. Because the focus of this work is on target selection and control, the
initial estimates were set to the truth values and the associated covariances
were set to be $\mathbf{P}_{0}=\beta\mathbf{I}_{n\times n}$ where $\beta$ is a
randomly generated parameter $10\leq\beta\leq 100$. If discrepancies in the
initial truth state and state estimate were considered, the complexity of the
problem increases dramatically. Elements of target search are introduced if
the discrepancy between truth position and estimated position grows
sufficiently large. This alone introduces the question: what bound on initial
estimate error is feasible for estimate convergence? This question along with
related questions are too broad in scope to be considered here and have
consequently been left for future work. However, the relevancy and
significance of tackling such a problem cannot be overstated.

Figure 2 depicts azimuth-elevation tracks of the deputy spacecraft and the
ability of the chief to track each deputy. The orientation of the chief
switches rather frequently at the beginning of the simulation while the
covariance matrices begin to settle below the specified threshold of
$\epsilon=0$ on the covariance entropy seen in Figure 3. To prevent
instantaneous switching, a switch is only triggered when the entropy of the
associated covariance matrix falls below this threshold, in and the entropy of
the covariance of any other deputy exceeds the same threshold. A time
hysteresis is also applied indicating the target agent cannot switch within
100 seconds of one another such that $\Delta=100$ from Algorithm 1. It is
shown that as the simulation progresses, the chief switches less frequently
because these threshold are crossed less often as the chief becomes more
certain of the trajectory of each deputy. Figure 4 is the resulting control
torque to track the reference trajectory. Uncharacteristically large spikes
are present at about 1,750 seconds when the azimuth of the reference
trajectory crosses from $\pi$ to $-\pi$, causing the chief’s attitude to wrap
around. This particular cases is present due to the parameterization of the
attitude, and for this reason future work is needed to adapt the presented
methodology to attitude parameterizations that do not present this
discontinuity such as rotation matrices and unit quaternions.

Figure 2: The MPC tracks the reference azimuth and elevation from the deputy
selected by the supervisory algorithm. Figure 3: Covariance entropy as the
chief observes each target deputy. It is shown that once each entropy value
falls below the threshold, the control algorithm is able to maintain the
entropy below the threshold. Figure 4: Applied torque for the 3 deputy case.
Figure 5: NMTs that generate the azimuth-elevation tracks.

The performance of the algorithm was also tested on a 10 deputy case, shown in
Figure 6, Figure 7, Figure 8. In Figure 6 the MPC is able to correctly track
the reference trajectory specified by the supervisor algorithm even with
increased deputy count. However, because of the increased number of deputies,
the chief is not able to linger on each deputy for as long and must quickly
move on to track the next deputy. Despite the increased workload, the chief is
still able to achieve the desired level of confidence in each state estimate
indicated by Figure 7. But, each entropy stays closer to the threshold when
compared to the 3 deputy case. Furthermore, the chief must apply a greater
amount of torque throughout the duration of the the mission seen in Figure 8.
This is a direct result of the supervisor algorithm greedily selecting the
next reference deputy based on the information entropy. It is possible to
include control effort in the supervisor algorithm to help mitigate the effort
applied by the chief.

Figure 6: The algorithms are capable of tracking the reference azimuth and
elevation track when more deputies require tracking. Figure 7: The algorithms
are also able to maintain the covariance entropy below the desired threshold,
even though it takes longer to achieve this level of estimator confidence.
Figure 8: Applied torque for the 10 deputy case.

##  6 Conclusions

This paper has provided a solution to the local catalog maintenance problem,
which involves a controllable chief spacecraft performing state estimation
across multiple deputy agents. The maintenance problem is, at its core, a
decision-making problem. Thus, this work proposes a hierarchical supervising
algorithm which autonomously decides which of the nearby deputies must be
observed. A model predictive controller is implemented to generate a
trajectory of torque inputs that are optimal with respect to minimizing the
deputies’ state uncertainties. We demonstrate through numerical simulation
that this supervisor scheme succeeds in handling both short-term and long-term
estimation requirements.

## References

  * [1] A. C. O. e. al., “Economic Benefits of the Global Positioning System (GPS), Final Report,” RTI Report Number 0215471., 2019.
  * [2] T. Sgobba and F. A. Allahdadi, “Chapter 8 - Orbital Operations Safety,” Safety Design for Space Operations (F. A. Allahdadi, I. Rongier, and P. D. Wilde, eds.), pp. 411–602, Oxford: Butterworth-Heinemann, 2013.
  * [3] K. L. Hobbs, J. B. Lyons, M. S. Feather, B. P. Bycroft, S. Phillips, M. Simon, M. Harter, K. Costello, Y. Gawdiak, and S. Paine, “Space Trusted Autonomy Readiness Levels,” IEEE Aerospace Conference, 2023, 10.48550/ARXIV.2210.09059.
  * [4] R. Opromolla, G. Fasano, G. Rufino, and M. Grassi, “Pose estimation for spacecraft relative navigation using model-based algorithms,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 53, No. 1, 2017, pp. 431–447.
  * [5] H. Patel, T. A. Lovell, S. Allgeier, R. Russell, and A. Sinclair, “Relative navigation for satellites in close proximity using angles-only observations,” AAS/AIAA Space Flight Mechanics Meeting, American Astronomical Soc. Washington, DC, 2012, pp. 12–202.
  * [6] D. C. Woffinden and D. K. Geller, “Observability Criteria for Angles-Only Navigation,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 45, No. 3, 2009, pp. 1194–1208, 10.1109/TAES.2009.5259193.
  * [7] D. K. Geller and I. Klein, “Angles-only navigation state observability during orbital proximity operations,” Journal of Guidance, Control, and Dynamics, Vol. 37, No. 6, 2014, pp. 1976–1983.
  * [8] Y. Hu, I. Sharf, and L. Chen, “Three-spacecraft autonomous orbit determination and observability analysis with inertial angles-only measurements,” Acta Astronautica, Vol. 170, 2020, pp. 106–121.
  * [9] M. L. Psiaki, “Autonomous Orbit Determination for Two Spacecraft from Relative Position Measurements,” Journal of Guidance, Control, and Dynamics, Vol. 22, No. 2, 1999, pp. 305–312.
  * [10] R. S. Erwin, P. Albuquerque, S. K. Jayaweera, and I. Hussein, “Dynamic sensor tasking for Space Situational Awareness,” Proceedings of the 2010 American Control Conference, 2010, pp. 1153–1158.
  * [11] B. Schutz, B. Tapley, and G. H. Born, Statistical orbit determination.  Elsevier, 2004.
  * [12] M. Sciotti, P. Besso, T. Flohrer, and H. Krag, “Low Earth Orbit objects tracking and orbit determination from ground-based phased array radar systems,” 2011 12th International Radar Symposium (IRS), IEEE, 2011, pp. 591–596.
  * [13] A. Herz and O. Logic, “SSA sensor tasking approach for improved orbit determination accuracies and more efficient use of ground assets,” 2017.
  * [14] A. D. Jaunzemis, M. J. Holzinger, and M. K. Jah, “Evidence-based sensor tasking for space domain awareness,” Advanced Maui Optical and Space Surveillance Technologies Conference, Maui Economic Development Board Maui, HI, 2016, p. 33.
  * [15] M. R. Endsley, “Toward a theory of situation awareness in dynamic systems,” Situational awareness, pp. 9–42, Routledge, 2017.
  * [16] M. Holzinger and M. Jah, “Challenges and potential in space domain awareness,” 2018.
  * [17] G. R. Frey, C. D. Petersen, F. A. Leve, I. V. Kolmanovsky, and A. R. Girard, “Constrained spacecraft relative motion planning exploiting periodic natural motion trajectories and invariance,” Journal of Guidance, Control, and Dynamics, Vol. 40, No. 12, 2017, pp. 3100–3115.
  * [18] J. L. Crassidis and J. L. Junkins, Optimal estimation of dynamic systems.  Chapman and Hall/CRC, 2004.
  * [19] C. W. Hays, T. Henderson, S. Phillips, and A. Soderlund, “Enabling Proliferated Space Sensor Awareness Constellations Through Topological Observability,” 2023 AIAA SciTech Forum, AIAA, 2023.
  * [20] A. Ly, M. Marsman, J. Verhagen, R. P. Grasman, and E.-J. Wagenmakers, “A tutorial on Fisher information,” Journal of Mathematical Psychology, Vol. 80, 2017, pp. 40–55.
  * [21] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi – A software framework for nonlinear optimization and optimal control,” Mathematical Programming Computation, Vol. 11, No. 1, 2019, pp. 1–36, 10.1007/s12532-018-0139-4.