ディープ・スペース・ネットワークのスケジューリング問題#
導入: ディープ・スペース・ネットワークとは？#
ディープ・スペース・ネットワーク（以下、DSN）とは、NASAが保有する電波望遠鏡のネットワークで、オーストラリア・スペイン・アメリカにまたがっています。
これらの電波望遠鏡は主に、宇宙に打ち上げられた宇宙探査機との交信、電波波長帯域での惑星観測ミッションに役立てられています。
しかし、宇宙探査機の増加やその交信の複雑化が進む昨今では、どのアンテナをどの宇宙探査機との交信に用いるか（あるいは、どのミッションに割り当てるか）のスケジューリングは困難を極めます。
近年、複数の電波望遠鏡をネットワーク化し、あるパターンに従って操作することで、電波信号の干渉を利用する電波干渉法が注目されています。
最も顕著な例は、M87銀河と天の川銀河の中心に存在する、超大質量ブラックホールの直接撮像に成功した、イベント・ホライズン・テレスコープ (EHT)です。
このような手法は今では主流となり、今後もスケジューリングの困難さは増していくことでしょう。
Guillaume et al. (2022)では、このスケジューリング問題をQUBO定式化し、D-Waveのハイブリッドソルバーを用いて解きました。今回はその定式化をJijModelingで実装し、JijZeptSolverで解いてみましょう。
DSNスケジューリング問題#
問題設定#
\(N\) 個のリクエストがあり、これを \(\mathcal{Q} = \{Q_1, Q_2, \dots, Q_N\}\) とします。あるリクエスト \(Q_n\) が \(M\) 個の(アンテナや装置などの地上の)リソースに割り当てることを考えます。このリソースの集合を \(\mathcal{S} = \{S_1, S_2, \dots, S_M\}\) とするとき、各地上のリソース \(S_m\) に対し、宇宙にある衛星が見える期間を \(\mathcal{V} = \{V_1, V_2, \dots, V_K\}\) とします。衛星が見える期間の開始時刻(rise time)と終了時刻(set time)をそれぞれrt, stと定義します。また、実際に衛星に信号を送信が開始できるようになるのはrtからわずかに遅れた時刻であり、送信ができるのはstからわずかに早い時刻までであるため、実際に送信が開始できる時刻と送信が終了する時刻をtn, tfと定義することとします。加えて、各リクエストに対して、要求された装置を準備する時間、装置を撤去するための時間が必要なので、それらをsu, tdと定義します。さらに、各リクエストにはトラッキングの継続時間drが定められていることとします。
これらをインスタンスデータとして、最適化を行います。DSNスケジュール最適化の目的は、全てのリクエストを満たすこと、言い換えればリクエストごとに正確に一つのアクションをスケジュールすることです。そのために、 \(x_{n, m, k, t}\) のようなバイナリ変数を用意しましょう。これは時刻 \(t\) からトラッキングを開始し、衛星が見える期間 \(k\) の間にリソース \(m\) を用いてリクエスト \(n\) を処理する場合に1、そうでない場合に0となるようなものです。
制約1: 全てのリクエストは決められた時間内に、どれかのアンテナで処理されなければならない#
制約1を数式で表現すると、以下のようになります。
実際にリクエスト送信を開始できるtnから、tf-drまでにリクエストの送信を終える必要があります。drだけトラッキング継続時間がかかるため、その分早くリクエストを送信しておかなければなりません。
制約2: プロジェクト間の衝突を避ける#
次の図に示すように、同じアンテナで同じ時間に2つのスケジュールが重なってはなりません。
上の図では、リクエスト \(i\) と \(j\)、 \(i\) と \(k\) は重なっていませんが、 \(j\) と \(k\) が重なってしまっています。このようなスケジュールの重なりは避けねばなりません。これを数式で表現すると以下のようになります。
JijModelingによる定式化#
以下では、JijModelingを用いて上記で定式化した数理モデルを実装します。
変数の定義#
式(1)および式(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 は各リクエストの各リソースと各「衛星が見える期間」についての処理開始可能時刻の範囲を表すタプルのリスト、 max_AT は処理開始可能時刻として考えなければならない時刻の最大値、 N はリクエストの総数、 M は地上のリソースの数、 K は宇宙にある衛星の見える期間の総数として定義しています。 su, dr, td はそれぞれ装置の準備時間、トラッキングの継続時間、装置の撤去時間を定義しています。そして、決定変数 \(x_{n, m, k, t}\) として x を定義しています。
制約#
次に、式(1)および式(2)を制約とする数理モデルを実装します。
@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]))
],
)
JijModelingでは、式(1)に含まれる \(\sum_{m, k, t}\) を jm.sum の内部でジェネレータ式を使って定義します。また、式(2)に含まれる複雑な条件式も、リスト内包表記と if 節を使って実装できます。
これで数理モデルの実装は完了です。正しく数理モデルが実装されているかをLaTeX表示を通して確認してみましょう。
problem
インスタンスの準備#
数理モデルの実装が完了したので、パラメーターに入力する値としてインスタンスデータを準備しましょう。今回の説明では以下のコードから生成されるインスタンスデータを利用します。
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 a array of transmission-on time and duration: [[(tn, tf - tn) for each viewperiod] 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 viewperiods
inst_M = inst_transmission_durations.shape[0]
inst_K = inst_transmission_durations.shape[1]
# compute a array of available time range [tn, tf - dr] (represented by
# (tn, (tf - dr) - tn + 1), where the second component is discretized time duration
# (end-inclusive) ralative to tn)
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}
JijZeptSolverで解く#
この問題をjijzept_solverで解いてみましょう。
import jijzept_solver
instance = problem.eval(instance_data)
solution = jijzept_solver.solve(instance, solve_limit_sec=1.0)
解の可視化#
最後に、得られた解を可視化してみましょう。
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()
上の図では、灰色が各リソースと宇宙探査機が通信できる時間を示しており、青色がトラッキング時間、黄色とマゼンタが各リクエストの準備時間と撤去時間を示しています。全ての青色（トラッキング時間）の区間が灰色（通信可能時間）に収まっているため、このスケジュールであればすべての通信を実現できることがわかります。