Deep Space Network Scheduling

Introduction: Deep Space Network (DSN)

NASA’s Deep Space Network (hereafter referred to as DSN) is a network of radio telescopes. These sites are Australia, Spain, and the United States. These radio telescopes are used primarily for communication with spacecraft launched into space and for planetary observation missions in the radio wavelength band. With the increasing number and complexity of spacecraft, communication scheduling with spacecraft is an extremely difficult problem. In recent years, radio interferometry, a technique that uses interference between signals by networking multiple radio telescopes and manipulating them according to a pattern, has attracted much attention. The most notable example is Event Horizon Telescope, which took the direct imaging of the supermassive black hole at the center of M87 and the Milky Way Galaxy.
Guillaume et al. (2022) formulated the DNS scheduling problem into QUBO and solved it using D-Wave’s hybrid solver. Here, we implement this problem using JijModeling and solve it with JijZeptSolver.

DSN scheduling problem

Problem setting

A schedule is generated after processing \(N\) requests \(\mathcal{Q} = \{Q_1, Q_2, \dots, Q_N\}\). A given request \(Q_n\) can utilize any of the \(M\) resources (antennas or equipments) prescribed for this request \(\mathcal{S} = \{S_1, S_2, \dots, S_M\}\). In turn, periods of visibility of a spacecraft from a particular ground station are calculated and define a set of \(K\) view periods \(\mathcal{V} = \{V_1, V_2, \dots, V_K \}\) for each resource \(S_m\). The times when a spacecraft becomes visible or invisible, respectively, the rise and set times, \(t^\mathrm{rt}\) and \(t^\mathrm{st}\), define the outer boundaries of a view period. The transmission to the spacecraft can begin (end) slightly later (earlier) than the rise and set time (respectively). These transmission on and off times, \(t^\mathrm{tn}\) and \(t^\mathrm{tf}\), delimit the inner boundaries within which a track can be scheduled. For each request, a setup and tear-down time period, \(\Delta t^\mathrm{su}\) and \(\Delta t^\mathrm{td}\), are given to prepare and remove (respectively) the equipment requested. Those times can occur anywhere between the rise and set times. In addition, a track duration \(\Delta t^\mathrm{dr}\) is prescribed for each request. An activity is defined as the period of time equal to the sum of the setup, track, and tear-down times, \(\Delta t^\mathrm{su}+\Delta t^\mathrm{dr}+\Delta t^\mathrm{td}\). It is the period of time during which an antenna can track a spacecraft and communicate with it. The goal of DSN schedule optimization is to process all requests; in other words, to schedule exactly one action per request. To this end, we use binary variables \(x_{n, m, k, t}\). These are equal to 1 when a track starts at time \(t\) within the view period \(k\) calculated for resource \(m\) of request \(n\), and \(0\) otherwise.

Constraint 1: Each request must be processed by one of the resources

\[ \sum_{m=1}^M \sum_{k=1}^K \sum_{t} x_{n, m, k, t} = 1 \qquad \forall n, t^\mathrm{tn} \leq t \leq t^\mathrm{tf}-\Delta t^\mathrm{dr} \tag{1} \]

Each request must be processed by one of the resources. The request must be completed by \(t^\mathrm{tf}- \Delta t^\mathrm{dr}\) from the time \(t^\mathrm{tn}\) actually starts sending the request. The tracking time \(\Delta t^\mathrm{dr}\) must be taken into account.

Constraint 2: Avoiding conflicts

In addition, the schedule should be devoid of conflicts as the figure below.

The figure shows requests \(i\) and \(j\), \(i\) and \(k\) do not conflict, but \(j\) and \(k\) overlap, indicating that they are in conflict. We have to avoid such conflicts. This can be expressed in the formula as follows:

\[\begin{split} \begin{align} & x_{n, m, k, t} x_{n', m', k', t'} = 0 \\ & \quad (Q_n \neq Q_{n'}, S_m = S_{m'}, t-\Delta t^\mathrm{su} \leq t'-\Delta t^\mathrm{su'} \leq t + \Delta t^\mathrm{dr} + \Delta t^\mathrm{td} \ \mathrm{or} \ t' - \Delta t^\mathrm{su'} \leq t-\Delta t^\mathrm{su} \leq t' + \Delta t^\mathrm{dr'} + \Delta t^\mathrm{td'}) \tag{2} \end{align} \end{split}\]

Implementation

Here, we implement a script that solves this problem using JijModeling and JijZeptSolver.

Define variables

We define the variables for eq.(1) and (2).

import jijmodeling as jm

problem = jm.Problem("Radio_telescope_scheduling")

AT = problem.Natural("AT", ndim=4, jagged=True)
max_AT = problem.DependentVar("max_AT", AT.max())
N = problem.Natural("N")
M = problem.Natural("M")
K = problem.Natural("K")
su = problem.Float("su", shape=(N,))
dr = problem.Float("dr", shape=(N,))
td = problem.Float("td", shape=(N,))
x = problem.BinaryVar("x", shape=(N, M, K, max_AT+1))

AT is the starting time of processing requests, max_AT is the maximum value of AT, N is the number of requests, M is the number of ground resources, and K is the number of viewperiods. su, dr, and td are the setup time \(\Delta t^\mathrm{su}\), tracking duration \(\Delta t^\mathrm{dr}\), and teardown time \(\Delta t^\mathrm{td}\) for each request, respectively. Also, we define binary variables x for \(x_{n, m, k, t}\).

@problem.update
def _(problem: jm.DecoratedProblem):
    problem += problem.Constraint(
        "onehot",
        [
            jm.sum(x[n1, m, k1, t1] for m in M for k1 in K for t1 in AT[n1, m, k1]) == 1
            for n1 in N
        ],
    )
    problem += problem.Constraint(
        "conflict",
        [
            x[n1, m, k1, t1] * x[n2, m, k2, t2] == 0
            for n1 in N
            for n2 in N
            if n2 != n1
            for m in M
            for k1 in K
            for k2 in K
            for t1 in AT[n1, m, k1]
            for t2 in AT[n2, m, k2]
            if (t1 - su[n1] <= t2 - su[n2]) & (t2 - su[n2] <= t1 + dr[n1] + td[n1])
            | ((t2 - su[n2] <= t1 - su[n1]) & (t1 - su[n1] <= t2 + dr[n2] + td[n2]))
        ],
    )

With JijModeling, we describe \(\sum_{m, k, t}\) in eq.(1) as jm.sum(x[n1, m, k1, t1] for m in M for k1 in K for t1 in AT[n1, m, k1]). To express all \(n_1, n_2\) combinations except \(n_1 = n_2\) in equation (2), we use if n2 != n1. Complex conditions like \(t_{n1} - \Delta t^\mathrm{su}_{n1} \leq t_{n2} - \Delta t^\mathrm{su}_{n2} \leq t_{n1} + \Delta t^\mathrm{dr}_{n1} + \Delta t^\mathrm{td}_{n1}\) are implemented with list comprehensions in Python.
With Jupyter Notebook, we can check the mathematical model implemented.

problem

\[\begin{split}\begin{array}{rl} \text{Problem}\colon &\text{Radio\_{}telescope\_{}scheduling}\\\displaystyle \min &\displaystyle 0\\&\\\text{s.t.}&\\&\begin{aligned} \text{conflict}&\quad \displaystyle {x}_{n1,m,k1,t1}\cdot {x}_{n2,m,k2,t2}=0\quad \forall \left(n1,n2,m,k1,k2,t1,t2\right)\;\text{s.t.}\;n1\in \left\{0,\ldots ,N-1\right\},n2\in \left\{0,\ldots ,N-1\right\},n2\neq n1,m\in \left\{0,\ldots ,M-1\right\},k1\in \left\{0,\ldots ,K-1\right\},k2\in \left\{0,\ldots ,K-1\right\},t1\in {AT}_{n1,m,k1},t2\in {AT}_{n2,m,k2},t1-{su}_{n1}\leq t2-{su}_{n2}\land t2-{su}_{n2}\leq t1+{dr}_{n1}+{td}_{n1}\lor t2-{su}_{n2}\leq t1-{su}_{n1}\land t1-{su}_{n1}\leq t2+{dr}_{n2}+{td}_{n2}\\\text{onehot}&\quad \displaystyle \sum _{m=0}^{M-1}{\sum _{k1=0}^{K-1}{\sum _{t1\in {AT}_{n1,m,k1}}{{x}_{n1,m,k1,t1}}}}=1\quad \forall n1\;\text{s.t.}\;n1\in \left\{0,\ldots ,N-1\right\}\end{aligned} \\&\\\text{where}&\\&\text{Decision Variables:}\\&\qquad \begin{alignedat}{2}x&\in \mathop{\mathrm{Array}}\left[N\times M\times K\times max\_{}AT+1;\left\{0, 1\right\}\right]&\quad &4\text{-dim binary variable}\\\end{alignedat}\\&\\&\text{Placeholders:}\\&\qquad \begin{alignedat}{2}AT&\in \mathop{\mathrm{JaggedArray}}\left[(-)\times (-)\times (-)\times (-);\mathbb{N}\right]&\quad &4\text{-dimensional array of placeholders with elements in }\mathbb{N}\\dr&\in \mathop{\mathrm{Array}}\left[N;\mathbb{R}\right]&\quad &1\text{-dimensional array of placeholders with elements in }\mathbb{R}\\K&\in \mathbb{N}&\quad &\text{A scalar placeholder in }\mathbb{N}\\M&\in \mathbb{N}&\quad &\text{A scalar placeholder in }\mathbb{N}\\N&\in \mathbb{N}&\quad &\text{A scalar placeholder in }\mathbb{N}\\su&\in \mathop{\mathrm{Array}}\left[N;\mathbb{R}\right]&\quad &1\text{-dimensional array of placeholders with elements in }\mathbb{R}\\td&\in \mathop{\mathrm{Array}}\left[N;\mathbb{R}\right]&\quad &1\text{-dimensional array of placeholders with elements in }\mathbb{R}\\\end{alignedat}\\&\\&\text{Dependent Variables:}\\&\qquad \begin{alignedat}{2}max\_{}AT&=\max _{\vec{\imath }}{{{\left(AT\right)}}_{\vec{\imath }}}&\quad &\in \mathbb{N}\\\end{alignedat}\end{array} \end{split}\]

Prepare an instance

Next, we create an instance.

import numpy as np

# set the number of requests
inst_N = 12
# set a list of set up time period: su
rng = np.random.default_rng(1234)
inst_su = rng.normal(2.0, 0.5, inst_N).tolist()
# set a list of track duration: dr
inst_dr = rng.normal(2.0, 0.5, inst_N).tolist()
# set a list of tear down time period: td
inst_td = rng.normal(1.5, 0.5, inst_N).tolist()
# set a array of transmission-on time: tn
inst_tn = np.array([[0, 6, 12, 18], [2, 8, 14, 20], [4, 10, 16, 22]])
# set a array of transmission-off time: tf
inst_tf = np.array([[4, 10, 16, 22], [6, 12, 18, 24], [8, 14, 20, 26]])
# get the number of resources and viewperiods
inst_M = inst_tn.shape[0]
inst_K = inst_tn.shape[1]
# compute a array of available time tn ≤ t ≤ tf - dr
inst_available = []
for n in range(inst_N):
    inst_available.append([])
    for i in range(inst_M):
        inst_available[n].append([])
        for j in range(inst_tf.shape[1]):
            inst_available[n][i].append(list(range(inst_tn[i, j], np.floor(inst_tf[i, j]-inst_dr[n]).astype(int))))
instance_data = {"AT": inst_available, "su": inst_su, "dr": inst_dr, "td": inst_td, "N": inst_N, "M": inst_M, "K": inst_K}

Solve by JijZeptSolver

We solve this problem using jijzept_solver.

import jijzept_solver

instance = problem.eval(instance_data, constraint_detection=False)
solution = jijzept_solver.solve(instance, solve_limit_sec=1.0)

Visualize the solution

Then, we visualize the solution obtained.

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

df = solution.decision_variables_df
x_indices = np.array(df[df["value"]==1]["subscripts"].to_list())
n_indices = x_indices[:, 0]
m_indices = x_indices[:, 1]
k_indices = x_indices[:, 2]
t_indices = x_indices[:, 3]
# make plot
fig, ax = plt.subplots()
# set x- and y-axis
ax.set_xlabel("Time")
ax.set_yticks(range(inst_M))
ax.set_yticklabels(["Resource {}".format(m) for m in range(inst_M)])
ax.get_xaxis().set_major_locator(ticker.MaxNLocator(integer=True))
# make bar plot for transmission using broken_barh
for n, m, k, t in zip(n_indices, m_indices, k_indices, t_indices):
    ax.broken_barh([(t, inst_dr[n])], (m-0.5, 1), color="dodgerblue")
    ax.broken_barh([(t-inst_su[n], inst_su[n])], (m-0.5, 1), color="gold")
    ax.broken_barh([(t+inst_dr[n], inst_td[n])], (m-0.5, 1), color="violet")
# make bar plot for available time
for m, (list_tn, list_tf) in enumerate(zip(inst_tn, inst_tf)):
    for tn, tf in zip(list_tn, list_tf):
        ax.broken_barh([(tn, tf-tn-1)], (m-0.5, 1), color="lightgray", alpha=0.4)
# make legend
ax.scatter([], [], color="lightgray", label="Transimission available", marker="s")
ax.scatter([], (), color="dodgerblue", label="Tracking", marker="s")
ax.scatter([], (), color="gold", label="Set up", marker="s")
ax.scatter([], (), color="violet", label="Tear down", marker="s")
ax.legend(bbox_to_anchor=(1.45, 1.0))
# show plot
plt.show() 

The gray regions indicate the times when each resource can communicate with the satellite. The blue bands represent the tracking. The yellow and magenta colors represent the setup time and teardown time for each request, respectively. We can see that all tracking durations overlap the gray communication available time.

References

• Guillaume et al., 2022, “Deep Space Network Scheduling Using Quantum Annealing”

• Event Horizon Telescope