整数長ジョブシーケンス問題#

こちらでは、Lucas, 2014, “Ising formulations of many NP problems”の 6.3. Job Sequencing with Integer Lengths を OpenJij と JijModeling、そしてJijModeling transpiler を用いて解く方法をご紹介します。

概要: 整数長ジョブシーケンス問題とは#

タスク1は実行するのに1時間、タスク2は実行に3時間、というように、整数の長さを持つタスクがいくつかあったとします。 これらを複数の実行するコンピュータに配分するとき、偏りを作ることなくコンピュータの実行時間を分散させるにはどのような組合せがあるでしょうか、というのを考える問題です。

具体例#

分かりやすくするために具体的に以下のような状況を考えてみましょう。

ここに10個のタスクと3個のコンピュータがあります。10個の仕事の長さはそれぞれ1,2,…,101, 2, \dots, 10とします。 これらのタスクをどのようにコンピュータに仕事を割り振れば仕事にかかる時間の最大値を最小化できるか考えます。 この場合、例えば1つ目のコンピュータには9, 10、2つ目には1, 2, 7, 8、3つ目には3, 4, 5, 6とするととなり、3つのコンピュータの実行時間の最大値は19となり、これが最適解です。

問題の一般化#

NN個のタスク{0,1,…,N−1}\{0, 1, \dots, N-1\}とMM個のコンピュータ{0,1,…,M−1}\{0, 1, \dots, M-1\}を考えましょう。各タスクの実行にかかる時間のリストをL={L0,L1,…,LN−1}\bm{L} = \{L_0, L_1, \dots, L_{N-1}\}とします。 jj番目のコンピュータで実行される仕事の集合をVjV_jとしたとき、コンピュータjjでタスクを終えるまでの時間はAj=∑i∈VjLiA_j = \sum_{i \in V_j} L_iとなります。 ii番目のタスクをコンピュータjjで行うことを表すバイナリ変数をxi,jx_{i, j}とします。

制約: タスクはどれか1つのコンピュータで実行されなければならない

例えば、タスク3をコンピュータ1と2の両方で実行することは許されません。これを数式にすると

(1)#∑j=0M−1xi,j=1(∀i∈{0,1,…,N−1}) \nonumber \sum_{j=0}^{M-1} x_{i, j} = 1 \quad (\forall i \in \{ 0, 1, \dots, N-1 \})

目的関数: コンピュータ1の実行時間と他の実行時間の差を小さくする

コンピュータ1の実行時間を基準とし、それと他のコンピュータの実行時間の差を最小にすることを考えます。これにより実行時間のばらつきが抑えられ、タスクが分散されるようになります。

(2)#min⁡{∑j=1M−1(A1−Aj)2} \nonumber \min\left\{ \sum_{j=1}^{M-1} (A_1 -A_j)^2\right\}

JijModelingによるモデル構築#

整数長ジョブシーケンス問題で用いる変数を定義#

式(1), (2)で用いられている変数を、以下のようにして定義しましょう。

import jijmodeling as jm


# defin variables
L = jm.Placeholder('L', dim=1)
N = L.shape[0]
M = jm.Placeholder('M')
x = jm.Binary('x', shape=(N, M))
i = jm.Element('i', (0, N))
j = jm.Element('j', (0, M))

L=jm.Placeholder('L', dim=1)でコンピュータに実行させるタスクの実行時間のリストを定義します。 そのリストの長さをN=L.shape[0]として定義します。Mはコンピュータの台数、xはバイナリ変数です。最後にxi,jx_{i, j}のように、変数の添字として使うものをi, jとして定義します。

制約の追加#

式(1)を制約として実装します。

# set problem
problem = jm.Problem('Integer Jobs')
# set constraint: job must be executed using a certain node
onehot = x[i, :]
problem += jm.Constraint('onehot', onehot==1, forall=i)

問題を作成し、そこに制約を追加しましょう。x[i, :]とすることでSum(j, x[i, j])を簡潔に実装することができます。

目的関数の追加#

式(2)の目的関数を実装しましょう。

# set objective function: minimize difference between node 0 and others
diffj = jm.Sum(i, L[i]*x[i, 0]) - jm.Sum(i, L[i]*x[i, j])
sumdiff2 = jm.Sum((j, j!=0), diffj*diffj)
problem += sumdiff2

diffjでA1−AjA_1 - A_jを計算し、それを2乗して総和を取ったものを制約とします。
実際に実装された数式をJupyter Notebookで表示してみましょう。

インスタンスの作成#

実際に実行するタスクなどを設定しましょう。

# set a list of jobs
inst_L = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
# set the number of Nodes
inst_M = 3
instance_data = {'L': inst_L, 'M': inst_M}

先程の具体例と同様に、{1,2,…,10}\{1, 2, \dots, 10\}の10個のタスクを3台のコンピュータで実行することを考えます。

未定乗数の設定#

整数長ジョブシーケンスには制約が一つあります。よってその制約の重みを設定する必要があります。 先程のConstraint部分で付けた名前と一致させるように、辞書型を用いて設定を行います。

# set multipliers
lam1 = 1.0
multipliers = {'onehot': lam1}

JijModeling transpilerによるPyQUBOへの変換#

ここまで行われてきた実装は、全てJijModelingによるものでした。 これをPyQUBOに変換することで、OpenJijはもちろん、他のソルバーを用いた組合せ最適化計算を行うことが可能になります。

from jijmodeling.transpiler.pyqubo import to_pyqubo

# convert to pyqubo
pyq_model, pyq_chache = to_pyqubo(problem, instance_data, {})
qubo, bias = pyq_model.compile().to_qubo(feed_dict=multipliers)

JijModelingで作成されたproblem、そして先ほど値を設定したinstance_dataを引数として、to_pyquboによりPyQUBOモデルを作成します。次にそれをコンパイルすることで、OpenJijなどで計算が可能なQUBOモデルにします。

OpenJijによる最適化計算の実行#

今回はOpenJijのシミュレーテッド・アニーリングを用いて、最適化問題を解くことにします。 それには以下のようにします。

import openjij as oj

# set sampler
sampler = oj.SASampler()
# solve problem
response = sampler.sample_qubo(qubo)

SASamplerを設定し、そのサンプラーに先程作成したQUBOモデルのquboを入力することで、計算結果が得られます。

デコードと解の表示#

返された計算結果をデコードし、解析を行いやすくします。

# decode solution
result = pyq_chache.decode(response)

このようにして得られた結果から、タスク実行が分散されている様子を見てみましょう。

# extract feasible solution
feasible = result.feasible()
# get the index of the lowest objective function
objectives = feasible.evaluation.objective
obs_dict = {i: objectives[i] for i in range(len(objectives))}
lowest = min(obs_dict, key=obs_dict.get)
# get indices of x = 1
indices, _, _ = feasible.record.solution['x'][lowest]
# get task number and execution node
tasks, nodes = indices
# get instance information
L = instance_data['L']
M = instance_data['M']
# initialize execution time
exec_time = [0] * M
# compute summation of execution time each nodes
for i, j in zip(tasks, nodes):
    exec_time[j] += L[i]
y_axis = range(0, max(exec_time)+1, 2)
node_names = [str(j) for j in range(M)]
fig = plt.figure()
plt.bar(node_names, exec_time)
plt.yticks(y_axis)
plt.xlabel('Computer numbers')
plt.ylabel('Execution time')
fig.savefig('integer_jobs.png')

すると以下のように、3つのコンピュータの実行時間がほぼ均等な解が得られます。