OpenJijで最適化問題を解く

OpenJijで最適化問題を解く#

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

3つ目のセクションでは、 ommx-openjij-adapterを使用します。ommx-open-jij-adapterはpipを使用してインストールできます。

pip install ommx-openjij-adapter

概要#

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出力する以外にはあまりできることはありません。インスタンスデータと組み合わせて、ソルバーの入力を生成することが必要になります。

まずは数理モデルの Placeholder に入力するインスタンスデータを用意し、 Interpreter オブジェクトに登録します。

Interpreter クラスのコンストラクタの引数に、以下のキーと値を持つ辞書を渡すことでインスタンスデータを登録できます：

  • キー：Placeholder オブジェクトの name プロパティに設定した文字列

  • 値：入力するデータ

instance_data = {
 "W": 100,                                  # ナップサックの耐荷重のデータ
 "v": [100, 90, 80, 70, 60, 50, 40, 30],    # アイテムの価値のデータ
 "w": [1, 5, 10, 20, 30, 40, 50, 60, 70]    # アイテムの重さのデータ
}
interpreter = jm.Interpreter(instance_data)

数理モデルをインスタンスに変換するには、Interpreter.eval_problem メソッドを使用します。インスタンスデータが登録された Interpreter オブジェクトの eval_problem メソッドに Problem オブジェクトを渡すと、その Problem オブジェクトが持つ Placeholder にインスタンスデータを入力され、インスタンスに変換されます:

instance = interpreter.eval_problem(problem)

Hint

Interpreter.eval_problem の返却値は ommx.v1.Instance オブジェクトです。OMMXについて詳しく知りたい場合はこちらを参照してください。

OpenJijサンプラーの入力値に変換できるOMMXインスタンスを取得できました。ommx-openjij-adapterで実際サンプリングを行ってみましょう。

from ommx_openjij_adapter import OMMXOpenJijSAAdapter

# Solve through SCIP and retrieve results as an ommx.v1.Solution
samples = OMMXOpenJijSAAdapter.sample(instance, num_reads=1)

print(f"目的関数の最適値:: {samples.best_feasible().objective}")

目的関数の最適値:: 337.75

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

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

次のステップ#