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 DSN 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#
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:
Modeling by JijModeling#
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")
N = problem.Natural("N")
M = problem.Natural("M")
K = problem.Natural("K")
# "Available Time ranges": Each tuple (s, d) in AT[n][m][k] represents the time interval [s, s+d) during which
# satellite tracking for the n-th request can be *commenced* by the m-th telescope for the k-th
# time slot (of the satellite associated with the n-th request).
AT = problem.Placeholder("AT", shape=(N, M, K), dtype=(jm.DataType.NATURAL, jm.DataType.NATURAL))
# The maximum feasible start time we need to consider
max_AT = problem.NamedExpr("max_AT", jm.map(lambda s, d: (s + d) - 1, AT).max())
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 a list of tuples representing the range of feasible start times for each request, resource, and satellite view period. max_AT is the maximum start time that must be considered, N is the number of requests, M is the number of ground resources, and K is the number of satellite view periods.
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, AT[n1, m, k1][0] + si]
for m in M
for k1 in K
for si in AT[n1, m, k1][1]
) == 1
for n1 in N
],
)
problem += problem.Constraint(
"conflict",
[
x[n1, m, k1, s1 + s1i] * x[n2, m, k2, s2 + s2i] == 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 s1, s1d in [AT[n1, m, k1]]
for s1i in s1d
for s2, s2d in [AT[n2, m, k2]]
for s2i in s2d
# Note: We'd like to write `for t1 in [s1 + s1i]` and `for t2 in [s2 + s2i]` here
# and replace `s1 + s1i` and `s2 + s2i` with `t1` and `t2`, respectively,
# but jijmodeling 2.5.0 crashes when we do so.
if (((s1 + s1i) - su[n1] <= (s2 + s2i) - su[n2])
& ((s2 + s2i) - su[n2] <= (s1 + s1i) + dr[n1] + td[n1]))
| (((s2 + s2i) - su[n2] <= (s1 + s1i) - su[n1])
& ((s1 + s1i) - su[n1] <= (s2 + s2i) + dr[n2] + td[n2]))
],
)
With JijModeling, we describe \(\sum_{m, k, t}\) in eq. (1) using a generator expression inside jm.sum.
The complex conditions in eq. (2) can also be implemented with list comprehensions and if clauses.
This completes the implementation of the mathematical model. Let’s check whether the model has been implemented correctly by displaying it in LaTeX.
problem
Prepare an instance#
Next, we create an instance.
import numpy as np
rng = np.random.default_rng(0)
# set the number of requests
inst_N = 12
# set a list of set up time period: su
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 an array of transmission-on times and durations: [[(tn, tf - tn) for each view period] for each resource]
inst_transmission_durations = np.array([
[(0, 3), (6, 3), (12, 3), (18, 3)],
[(2, 3), (8, 3), (14, 3), (20, 3)],
[(4, 3), (10, 3), (16, 3), (22, 3)]]
)
# get the number of resources and view periods
inst_M = inst_transmission_durations.shape[0]
inst_K = inst_transmission_durations.shape[1]
# compute an array of available time ranges [tn, tf - dr] represented by
# (tn, (tf - dr) - tn + 1), where the second component is the discretized duration
# relative to tn and includes the end point
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_K):
tx_start_lower_bound, tx_duration = inst_transmission_durations[i, j]
tx_start_range_width = int(tx_duration - inst_dr[n]) + 1 # +1 for end-inclusiveness
inst_available[n][i].append((int(tx_start_lower_bound), tx_start_range_width))
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)
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, td_m in enumerate(inst_transmission_durations):
for tn, tf_minus_tn in td_m:
ax.broken_barh([(tn, tf_minus_tn)], (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 blue tracking intervals fit within the gray communication-available intervals.