ディープ・スペース・ネットワークのスケジューリング問題

導入: ディープ・ペース・ネットワークとは？

ディープ・スペース・ネットワーク（以下、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を数式で表現すると、以下のようになります。

\[ \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} \]

実際にリクエスト送信を開始できるtnから、tf-drまでにリクエストの送信を終える必要があります。drだけトラッキング継続時間がかかるため、その分早くリクエストを送信しておかなければなりません。

制約2: プロジェクト間の衝突を避ける

次の図に示すように、同じアンテナで同じ時間に2つのスケジュールが重なってはなりません。

上の図では、リクエスト \(i\) と \(j\)、 \(i\) と \(k\) は重なっていませんが、 \(j\) と \(k\) が重なってしまっています。このようなスケジュールの重なりは避けねばなりません。これを数式で表現すると以下のようになります。

\[\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}\]

実装しましょう

以下では、JijModelingを用いて上記で定式化した数理モデルを実装します。

変数の定義

式(1)および式(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 はリクエストの処理開始時間、 max_AT は 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, 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]))
        ],
    )

JijModelingでは、式(1)に含まれる \(\sum_{m, k, t}\) を jm.sum(M, lambda m: ...) のように定義します。また、式(2)に含まれる複雑な条件式も、lambda関数とjm.filterを使って実装することができます。

これで数理モデルの実装は完了です。正しく数理モデルが実装されているかをLaTeX表示を通して確認してみましょう。

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}\]

インスタンスの準備

数理モデルの実装が完了したので、パラメーターに入力する値としてインスタンスデータを準備しましょう。今回の説明では以下のコードから生成されるインスタンスデータを利用します。

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}

JijZeptSolverで解く

この問題をjijzept_solverで解いてみましょう。

import jijzept_solver

instance = problem.eval(instance_data, constraint_detection=False)
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, (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() 

上の図では、灰色が各リソースと宇宙探査機が通信できる時間を示しており、青色がトラッキング時間、黄色とマゼンタが各リクエストの準備時間と撤去時間を示しています。灰色（通信時間）が青色（トラッキング時間）にすべて塗りつぶされているため、このスケジュールであればすべての通信を実現できることがわかります。

参考文献

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

• イベント・ホライズン・テレスコープ