Cislunar Satellite Constellation Design via Integer Linear Programming

Suppose there is target point in space between the Earth and Moon, denoted by $\boldsymbol{p}$ and an observer at position $\boldsymbol{r}$. The target point’s observability can be quantified by an accessibility measure, $g$. A simple measure could be the distance between observer and target, $g(\boldsymbol{p},\boldsymbol{r})=\|\boldsymbol{p}-\boldsymbol{r}\|$. A more sophisticated measure is the apparent magnitude of the object, which is a function of the position of the sun as well, $g=g(\boldsymbol{p},\boldsymbol{r},\boldsymbol{s})$ [14].

Suppose the observer is on an orbit with period $T$, time is discretized for computational ease, and $\boldsymbol{p}$ is a static target point in space. At time step $n$, $\boldsymbol{r}_{n}$ is the observer’s position, $\boldsymbol{s}_{n}$ the sun’s position, and $a[n]=g(\boldsymbol{p},\boldsymbol{r}_{n},\boldsymbol{s}_{n})$ is the accessibility measure. If the orbit is discretized into $L$ equal time steps, then $a$ is a $L\times 1$ vector containing the accessibility measures along the observer’s orbit. In practice, the target point will be in one of two states, [observable] or [not observable]. In other words, the continuous accessibility measure must be transformed into a binary one. To achieve this, we introduce the accessibility profile, $v$, distinct from the accessibility measure vector, $a$. Define a threshold value, $M$, for the accessibility measure. Then,

At time step $n$, $v[n]$ indicates whether the observer can view the target point.

Instead of a single observer, suppose we have a set (called a constellation) of observers on the same orbit of period $T$ again discretized into L time steps. The constellation pattern vector $x$ captures the phasing of the observers in the constellation relative to a seed satellite. At time step $n=0$, if the seed satellite is at position $\boldsymbol{r}$, then after $n_{k}$ time steps, the $k$-th observer will be at position $\boldsymbol{r}$. The continuous-time time delay for the $k$-th observer is $\Delta t_{k}=n_{k}*\frac{T}{L}$. This is illustrated below in Figure 1.

The accessibility profile is a vector of binary variables describing when the target is in view of the observer as the observer travels along its orbit. The coverage timeline, denoted by $b$, is a $L\times 1$ vector that extends the accessibility profile concept to multiple observers. At time step $n$, $b[n]$ describes the number of satellites with view access to the target point. Instead of a vector of binary variables, the coverage timeline is a vector of non-negative integers.

It turns out that a circular convolution between the constellation pattern vector $x$ and the accessibility profile $v$ returns the coverage timeline $b$. The circular nature of the convolution is a direct consequence of the periodicity of the observers’ orbit. As such, $b$ can be computed from $x$ and $v$:

Alternatively, the circular convolution can be conveniently expressed as a matrix multiplication between an $L\times L$ circulant matrix, $V$, and the constellation pattern vector, $x$. Here, the $i$-th column of $V$ is produced by circularly shifting $v$ by $i$ elements.

Expanded in matrix form and using 1-based indexing,

Often times, we are interested in satisfying a coverage requirement, $f$, which demands a different number of observers at different time steps. The goal, then, is to place a set of observers on the orbit (i.e. find a constellation pattern $x$) such that the resulting coverage timeline $b$ meets or exceeds the coverage requirement $f$. Stated mathematically,

In general, there can be many constellation pattern vectors that satisfy the coverage requirement. However, it is often beneficial to minimize the total number of observers in the constellation while still satisfying the coverage requirement, since fewer observers reduce launch and maintenance costs.

The problem above of finding a constellation pattern $x$ is now a problem of finding the optimal constellation pattern $x^{*}$ that minimizes the number of observers in the constellation. The following integer linear program presents this mathematically.

The objective is to minimize the total number of satellites, hence the objective $\boldsymbol{1}^{T}x$, which sums all the entries in $x$.

The above formulation is linear, and thus allows extensions to multiple target points and multiple constellations.

As traffic between the Earth and the Moon is anticipated to increase, an attractive choice is a transfer orbit from L1 to geosynchronous Earth orbit (GEO). Figure 3 illustrates this transfer. It consists of a manifold insertion from L1 followed by a high impulse manuever into a connecting orbit to GEO. The transfer takes $\approx$ 23.12 days. The discretized position history (i.e. the blue circles in Figure 3) are target points comprising the set $\mathcal{J}$. Note that although we define $\mathcal{J}$ as a set of points, it is actually a sequence. Target point 1 is followed by target point 2 and so on, tracking the position of a hypothetical target on this path from L1 to GEO.

There are a variety of periodic orbits in the CR3BP. Of these, the most attractive options offer close approaches to L1 and GEO while moving through the majority of cislunar space. Furthermore, as per an AFRL objective to monitor a 30 degree cone emanating from the Earth, orbits in this region are also attractive. We choose a set of 6 orbits inspired by Gupta et al.[12], which are plotted in Figure 4 below. The initial conditions are shown in Table 3 in the appendix. All orbits were generated via a single-shooting differential correction algorithm with initial guesses sampled from data provided by NASA JPL[15] and Broucke. [16]. All orbits have period $T$ = 6.45 TU $\approx$ 28.008 days. It should be noted that the L1 Lyapunov and L2 Halo have a fundamental period of $\frac{T}{2}$ = 3.225 TU which means they are also $T$-periodic. The circular restricted problem does not admit many orbits with period equal to the synodic period of the Earth-Moon system (approximately 29.05 days). As a consequence, while the orbits themselves repeat, the illumination conditions will not. As such, the results of the formalism apply on a case-by-case basis depending on the sun’s initial phase (in radians) relative to the Earth-Moon axis, denoted as $\phi_{0}$. We stress that the candidate orbits are chosen to have periods as close as possible to the synodic period and are still resonant in the CR3BP sense. This ensures the difference in illumination conditions is small after $t=6.45$ TU compared to a randomly chosen constellation period. We aimed to approximately minimize this difference while preserving candidate orbit diversity. In this work we choose $\phi_{0}$ = 0 for optimization (Earth-Moon-Sun are collinear at simulation start). However, we investigate the robustness of the constellation to different initial phase angles later in the Results section.

To quantify the the observability of a target from an observer’s position, we utilize the target point’s apparent magnitude.[14]^,[9]Consider an observer at position $\boldsymbol{r}$, a spherical target with diameter $d$ at position $\boldsymbol{p}$, and the sun at position $\boldsymbol{s}$. Define the solar phase angle $\psi$,

It should be noted that $\psi$ is the familiar angle used in the dot product formula between two vectors, taking on values between 0 and $\pi$. The diffuse phase angle function can be written as a function of $\psi$,

On the interval (0, $\pi$), $p_{\text{diff}}$ is a strictly decreasing function. The apparent magnitude of a target is given by,

Where $\zeta=\|\boldsymbol{p}-\boldsymbol{r}\|$ is the distance between the oberver and target, and $a_{\text{spec}},a_{\text{diff}}$ are the specular and diffuse reflection coefficients of said target, respectively. Values can be found in Table 2 in the appendix. It should be noted that the solar phase angle $\psi$ is distinct from the initial phase angle $\phi_{0}$ described in the Observer Orbits section above. $\phi_{0}$ descibes the initial orientation of the Sun relative to the Earth-Moon axis. $\psi$ describes the orientation of an observer, a target point, and the Sun. The value of $g$ increases with $\psi$, meaning solar phase angles closer to $\pi$ result in an unfavorably dim target. Every unit increase in apparent magnitude corresponds to a 2.5 times dimmer target. In addition to this, we consider the possibility that the Earth and Moon may occlude the target. If this is the case, the accessibility measure is set to $\infty$, indicating an invisible target point. In practice, the measure will be set to the largest floating point number (often denoted REALMAX) available on the machine. We choose a threshold value of $M=17$ as our access measure cutoff. Together, $g$ and $M$ define our accessibility metric and threshold, respectively.

We seek a set $\mathcal{F}=\{f_{1},f_{2},\cdots,f_{\lvert\mathcal{J}\rvert}\}$ of coverage requirements, one for each target point in the set $\mathcal{J}$. As our object of interest travels along the transfer, it travels through each target point in order. Thus, the coverage requirement for target point $j+1$ should be that of target point $j$, but shifted (in fact, circularly shifted due to the periodicity of observer orbits) forward by one index. Mathematically,

Where $[\ ]$ symbolizes 0-based indexing into the coverage requirement array $f_{j}$. This recursive relation is convenient because we are only required to define the coverage requirement $f_{1}$, with subsequent $f_{j}$ defined by Eqn. 1. An example set $\mathcal{F}=\{f_{1},f_{2},f_{3}\}$ is illustrated below in Figure 5.

The constellation shows impressive robustness to delays in departure. Figures 9 and 9 show binary indicator maps for constellations optimized for $N=1$ and $N=16$ departure windows. Figure 15 in the appendix shows the maps for each $N\in\{1,2,4,8,16\}$ requested departure windows. Areas shaded in green indicate a pair (Time Step, Target Point) which is observable by at least one satellite in the constellation. Areas shaded in black are not observable. The red lines (shown as diagonal streaks across Figures 9 and 9) illustrate the spatio-temporal coverage requirement sent to the optimizer. These regions are always covered, as required by the formulation. We see that the majority of the optimizer’s work attempts to cover the target points closer to GEO (target point No. 300 and onwards). As seen in Figure 3, these points are spaced further apart (the target is moving quickly in this region), naturally making it more difficult to observe. With increasing $N$, we make progress in closing the gap, with most of the binary map displaying green for $N=16$.

We acknowledged earlier that because our candidate orbits are not resonant with the Earth-Moon synodic period, after one simulation period, the illumination conditions will differ slightly at the start of the second epoch ($t=6.45$ TU). The orbits repeat but now with a slightly different initial sun phase angle, $\phi_{0}$. As discussed above, we chose an initial sun phase angle of 0 rad, meaning that at the start of the simulation, the Sun is on the positive x-axis at distance 1 AU. The constellations are optimized for this phase angle, but they show impressive robustness against various initial sun phase angles. We investigated numerous $\phi_{0}$ between 0 and 2$\pi$, given a constellation already optimized for $N=2$ departure windows (see Figure 15 in the appendix for an illustration of the constellation). Figure 11 shows that despite changes in $\phi_{0}$, the constellation still satisfies a majority of the coverage requirement. As such, the constellation’s effectiveness extends beyond just one simulation period. For example, if mission requirements mandate a 70$\%$ coverage requirement satisfaction, then our constellation meets the mandate.

Furthermore, Figure 11 shows that there exist several values of $\phi_{0}$ still satisfying $>$ 90$\%$ of the coverage requirement. Figure 11 shows the coverage map for one such angle, $\phi_{0}=162^{\circ}$. Not only is 97$\%$ of the coverage requirement satisfied, the impressive ratio of green to black shows that robustness to delay is still present.

Of all the investigated $\phi_{0}$ in the range [0, 2$\pi$], $\phi_{0}=250^{\circ}$ was the worst performing constellation. However, it still satisfied 71.48 $\%$ of the coverage requirement. Figure 13 shows the coverage map for this choice of $\phi_{0}$, illustrating that despite only satisfying 71.48 $\%$ of the coverage requirement, the constellation still boasts impressive robustness to delays (shown via the ratio of green to black). There is a band of invisibility (appearing as an hourglass figure in black) in the approximate time step range 70 to 130. This is due to unfavorable solar phase angles $\psi$ that are incurred in this region. Figure 13 shows the position of the observers as well as the direction of the Sun’s rays at time step 100. We see that the positions of the Sun, target points, and observer on the Lyapunov orbit result in solar phase angle $>\frac{\pi}{2}$. In this scenario the Sun’s rays are too direct, resulting a poorly illuminated target, making it difficult for optical sensors to monitor.