OpenJijで最適化問題を解く

OpenJijで最適化問題を解く#

jijmodelingの使い方を理解するために基本的な例を見ていきましょう。以下では、シンプルな数理モデルを作成し、それを変換し、ソルバーで実行する手順を説明します。最初の2つのセクションでは、jijmodelingだけで十分ですが、LaTeX出力を簡単に確認するためにJupyterノートブックを使用することをお勧めします。

3つ目のセクションでは、JijModelingTranspilerとOpenJijを使用します。JijModelingTranspilerとOpenJijはpipを使用してインストールできます。

pip install jijmodeling-transpiler openjij

概要#

jijmodelingを記述方法は数理最適化に精通している人にとって自然に感じられるはずです。

jijmodelingのクラス（BinaryVarやPlaceholderなど）を基本的な演算を使用して組み合わせることで数式が表現できます。___Varクラスはさまざまな種類の決定変数を指します。Placeholderは後で指定される定数や値を表します。つまり、問題を抽象化し、インスタンスデータとしてマークしたいものを表しています。もちろん、jijmodelingによる数理モデルの構築には数値リテラルも使用できます。

import jijmodeling as jm

x = jm.BinaryVar("x")
y = jm.IntegerVar("y", lower_bound=1, upper_bound=10)
n = jm.Placeholder("n")
exp = x * (y ** 2) * n

上記のプレースホルダーと変数は0次元（スカラー）ですが、これらのクラスを使用して配列や多次元変数および定数を表現することもできます（後で詳しく説明します）。

jijmodelingでは、上記のexpのような式をProblemオブジェクトに加えることで数理モデルを構築し、数理モデル全体を表現します。制約は、比較式をラップするConstraintクラスによって定義されます（Constraintでは<=、==、および>=のみがサポートされています）。

a = jm.IntegerVar("a", lower_bound=5, upper_bound=20)
b = jm.IntegerVar("b", lower_bound=1, upper_bound=20)
c = jm.IntegerVar("c", lower_bound=20, upper_bound=30)
n = jm.Placeholder("n")

problem = jm.Problem("my problem")
problem += n * (a + b + c)
problem += jm.Constraint("c1", 2 * (b + c) <= 75)
problem += jm.Constraint("c2", a + b <= 40)
problem
\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{my problem} & & \\& & \min \quad \displaystyle n \cdot \left(a + b + c\right) & \\\text{{s.t.}} & & & \\ & \text{c1} & \displaystyle 2 \cdot \left(b + c\right) \leq 75 & \\ & \text{c2} & \displaystyle a + b \leq 40 & \\\text{{where}} & & & \\& a & 0\text{-dim integer variable}\\ & & \text{lower bound: }5 & \\ & & \text{upper bound: }20 & \\& b & 0\text{-dim integer variable}\\ & & \text{lower bound: }1 & \\ & & \text{upper bound: }20 & \\& c & 0\text{-dim integer variable}\\ & & \text{lower bound: }20 & \\ & & \text{upper bound: }30 & \\\end{array}\end{split}\]

モデルの作成#

一般的なナップサック問題をどのようにモデル化するか見てみましょう。

ナップサック問題とは、それぞれ価値と重量が設定された\(N\)個のアイテムがあるとき、重量制限\(W\)の範囲内で持ち帰れるアイテムの価値を最大化する問題です。以下では、アイテム\(i\)を持ち帰るかどうかをバイナリ変数\(x_{i}\)で表現し、アイテム\(i\)の重量を\(w_i\)、アイテム\(i\)の価値を\(v_i\)で表現するものとします。

まず、これらの値をjijmodelingで定義してみましょう。

# W: 重量制限（持ち帰れる重量の上限値）
W = jm.Placeholder("W")
# v_i: アイテムiの価値
values = jm.Placeholder("v", ndim=1) 
# w_i: アイテムiの重量
weights = jm.Placeholder("w", ndim=1) 
# Nはプレースホルダーの0次元方向のサイズに基づいて自動的に決まるため、
# 以下のように`len_at`メソッドを使って定義することができる。
# また、LaTeX出力のために、ここでは`latex`パラメータを指定しておく。
N = values.len_at(0, latex="N")

# x_i: アイテムiを持ち帰るかどうかを表すバイナリ変数
# ここでは、`shape`パラメータを使って、N個の要素を持つベクトルを定義している。
x = jm.BinaryVar("x", shape=(N,))

# 総和を取るためにインデックスを定義しておく。
i = jm.Element("i", (0, N))

上記のコードの最後に定義したElementを利用してjijmodeling.sumを使用すると、シグマ記法に似たスタイルで総和を記述できます。上記のコードでは、インデックスiを0を含み\(N\)を含まないの区間として定義しています。これを事前に記述するのは奇妙に感じるかもしれませんが、そうすることで再利用が可能になり利便性が向上します。

ナップサック問題では、重量制限を満たしながら持ち帰るアイテムの価値を最大化する必要があります。これをjijmodelingで表現すると以下のようになります。

problem = jm.Problem("knapsack", sense = jm.ProblemSense.MAXIMIZE)
problem += jm.sum(i, values[i] * x[i])
problem += jm.Constraint("weight limit", jm.sum(i, weights[i] * x[i]) <= W)
problem
\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{knapsack} & & \\& & \max \quad \displaystyle \sum_{i = 0}^{N - 1} v_{i} \cdot x_{i} & \\\text{{s.t.}} & & & \\ & \text{weight limit} & \displaystyle \sum_{i = 0}^{N - 1} w_{i} \cdot x_{i} \leq W & \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}\\\end{array}\end{split}\]

jijmodelingの式は通常のPythonオブジェクトのように変数に格納できます。これにより、大規模な問題に取り組む際には複雑な式を小さなものから構築することができ、後で理解・修正しやすくなります。ナップサック問題のような小さい問題ではあまり役に立ちませんが、例として、上記のコードを次のように書き直すことができます。

chosen_v = values[i] * x[i]
chosen_w = weights[i] * x[i]
sum_of_values = jm.sum(i, chosen_v)
sum_of_weights = jm.sum(i, chosen_w)
weight_below_limit = sum_of_weights <= W

problem = jm.Problem("knapsack", sense = jm.ProblemSense.MAXIMIZE)
problem += sum_of_values
problem += jm.Constraint("weight limit", weight_below_limit)
problem
\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{knapsack} & & \\& & \max \quad \displaystyle \sum_{i = 0}^{N - 1} v_{i} \cdot x_{i} & \\\text{{s.t.}} & & & \\ & \text{weight limit} & \displaystyle \sum_{i = 0}^{N - 1} w_{i} \cdot x_{i} \leq W & \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}\\\end{array}\end{split}\]

これらの2つのコードが表すモデルは同等です。あなたの好みや利便に合わせてモデルを記述してみてください。

モデルを解く#

前述のコードによりモデルを構築することができましたが、jijmodelingだけではLaTeX出力する以外にはあまりできることはありません。インスタンスデータと組み合わせて、ソルバーの入力を生成すること必要になります。例として、無料で利用できるJijModelingTranspilerを使用して数理モデルを変換し、OpenJijを使用してモデルを解く方法を見てみましょう。なお、JijZeptで提供されているツールもjijmodelingによって構築された数理モデルを入力として受け入れるようになっています。

本説は簡易的なデモンストレーションを目的としているため、JijModelingTranspilerやOpenJijの説明は行いません。詳細については、それぞれのドキュメントを参照してください。

以下では、JijModelingTranspilerを使用して、jijmodelingで構築した数理モデルとインスタンスデータを渡してQUBOに変換することとします。

from jijmodeling_transpiler.core import compile_model
from jijmodeling_transpiler.core.pubo import transpile_to_pubo

data = {
 "W": 100,
 "v": [100, 90, 80, 70, 60, 50, 40, 30],
 "w": [1, 5, 10, 20, 30, 40, 50, 60, 70]
}

compiled = compile_model(problem, data)
qubo, _ = transpile_to_pubo(compiled).get_qubo_dict()

これでOpenJijが入力として受け入れることができるQUBOを取得できました。QUBOを手作業で書き出すことはミスが非常に発生しやすいのですが、jijmodelingおよびjijmodeling_transpilerを利用すれば、そのようなミスを無くすことができ、検証が容易になります。では、実際にQUBOをopenjijに入力して問題を解いてみましょう。

import openjij as oj

sampler = oj.SASampler()
response = sampler.sample_qubo(qubo, num_reads=1)
response.first

Sample(sample={0: 1, 1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 0, 7: 0}, energy=-5.501111111111113, num_occurrences=1)

上記のコードはopenjijのシミュレーテッドアニーリングを使用しており、num_reads=1は1回だけサンプリングすることを示しています。num_readsの値を増やすことで複数回サンプリングし、レスポンスオブジェクトを使用してさまざまな結果を探索することができます。しかし、今回の問題では、すべてのサンプルが最適解に到達するため、ここでは1回のサンプルを行い、見つかった「最良」の解を見てみることとします。

Sample(sample={0: 1, 1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 0, 7: 0}, energy=-5.501111111111113, num_occurrences=1)

上記のレスポンスオブジェクトは各決定変数に対してopenjijが求めた値を示しています。JijModelingTranspilerやOpenJijを使用して結果をより良く処理および視覚化したり、同じ数理モデルを異なる目的で再利用したりするために、さらに多くのことができますが、それについてはそれぞれのドキュメントページを参照してください。

次のステップ#