JijModelingの発展的な使い方#
概要#
このドキュメントでは、JijZept Toolsを使ってJijModelingを発展的に使う方法について説明します。
JijModelingは、多くの数理最適化モデラーと異なり、入力された数式の構造をそのまま保持しています。これにより、数式の構造を自在に変形させることができ、この特性を応用した手法がJijZept Toolsに実装されています。
以下では、JijZept Toolsを用いた手法を紹介します。
本ドキュメントは、JijModelingの発展的な使い方に言及しています。JijModelingの基本的な使い方については、こちら を参照してください。
実行可能解が存在しない問題の、常に実行可能な問題への変換#
数理最適化において、対象となる問題に実行可能解が存在しない場合、どの制約が悪かったのかを特定することは一般に難しいです。
このような状況において、slack変数を導入することで、実行可能解が存在しない問題を常に実行可能な問題に変換することができ、さらにslack変数の値を見ることにより、どの制約が悪さをしているのかの見当をつけることができます。
実際に具体例を見てみましょう。
以下のJijModelingにより定式化された問題を考えます。
import jijmodeling as jm
# 問題の初期化
problem = jm.Problem('Binary_LP')
# パラメータの定義
A = jm.Placeholder('A', ndim=2, description='制約行列')
b = jm.Placeholder('b', ndim=1, description='制約の右辺ベクトル')
c = jm.Placeholder('c', ndim=1, description='コスト係数ベクトル')
# インデックスの定義j
m = A.len_at(0)
m.set_latex("M") # 行数
i = jm.Element('i', belong_to=(0, m), description='制約の行インデックス')
n = A.len_at(1)
n.set_latex("N") # 列数
j = jm.Element('j', belong_to=(0, n), description='変数の列インデックス')
# バイナリ変数の定義
x = jm.BinaryVar('x', shape=(n,), description='選択/非選択を表すバイナリ変数')
# 目的関数の定義
problem += jm.sum(j, c[j] * x[j]) # c^T x の最小化
# 制約の定義
problem += jm.Constraint(
'linear_constraint',
jm.sum(j, A[i, j] * x[j]) >= b[i],
forall=i
)
problem
\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{Binary\_LP} & & \\& & \min \quad \displaystyle \sum_{j = 0}^{N - 1} c_{j} \cdot x_{j} & \\\text{{s.t.}} & & & \\ & \text{linear\_constraint} & \displaystyle \sum_{j = 0}^{N - 1} A_{i, j} \cdot x_{j} \geq b_{i} & \forall i \in \left\{0,\ldots,M - 1\right\} \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}& \text{選択/非選択を表すバイナリ変数}\\\end{array}\end{split}\]
インスタンスデータを以下のように定義し、実際に求解してみます。
import numpy as np
from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter
instance_data = {
'A': np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
'b': np.array([10, 1, 20]),
'c': np.array([7, 8, 9])
}
interpreter = jm.Interpreter(instance_data)
ommx_instance = interpreter.eval_problem(problem)
# OMMXPySCIPOptAdapterを使用して問題を解く
solution = OMMXPySCIPOptAdapter.solve(ommx_instance)
/home/runner/work/JijZeptTools/JijZeptTools/.venv/lib/python3.10/site-packages/ommx_pyscipopt_adapter/adapter.py:30: UserWarning: linked SCIP 9.02 is not recommended for this version of PySCIPOpt - use version 9.2.1
self.model = pyscipopt.Model()
---------------------------------------------------------------------------
InfeasibleDetected Traceback (most recent call last)
Cell In[2], line 13
10 ommx_instance = interpreter.eval_problem(problem)
12 # OMMXPySCIPOptAdapterを使用して問題を解く
---> 13 solution = OMMXPySCIPOptAdapter.solve(ommx_instance)
File ~/work/JijZeptTools/JijZeptTools/.venv/lib/python3.10/site-packages/ommx_pyscipopt_adapter/adapter.py:132, in OMMXPySCIPOptAdapter.solve(ommx_instance, use_sos1)
130 model = adapter.solver_input
131 model.optimize()
--> 132 return adapter.decode(model)
File ~/work/JijZeptTools/JijZeptTools/.venv/lib/python3.10/site-packages/ommx_pyscipopt_adapter/adapter.py:182, in OMMXPySCIPOptAdapter.decode(self, data)
180 # there appears to be no enum for this in pyscipopt
181 if data.getStatus() == "infeasible":
--> 182 raise InfeasibleDetected("Model was infeasible")
184 if data.getStatus() == "unbounded":
185 raise UnboundedDetected("Model was unbounded")
InfeasibleDetected: Model was infeasible
実行不可能のエラーが出てしまいました。このままではどの制約が悪いのか分かりません。
ここで、JijZept Toolsを使って、slack変数を導入してみましょう。
以下の1行のコードを実行すれば良いです。
from jijzepttools.modeling.replace.problem import generate_always_feasible_problem_forall
replaced_problem = generate_always_feasible_problem_forall(problem)
replaced_problem
\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{Binary\_LP} & & \\& & \min \quad \displaystyle \sum_{j = 0}^{N - 1} c_{j} \cdot x_{j} + bigM \cdot \sum_{i = 0}^{M - 1} \left(v0_{i} + w0_{i}\right) & \\\text{{s.t.}} & & & \\ & \text{linear\_constraint} & \displaystyle \sum_{j = 0}^{N - 1} A_{i, j} \cdot x_{j} + v0_{i} - w0_{i} \geq b_{i} & \forall i \in \left\{0,\ldots,M - 1\right\} \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}& \text{選択/非選択を表すバイナリ変数}\\& v0 & 1\text{-dim integer variable}\\ & & \text{lower bound: }0 & \\ & & \text{upper bound: }100000000 & \\& w0 & 1\text{-dim integer variable}\\ & & \text{lower bound: }0 & \\ & & \text{upper bound: }100000000 & \\\end{array}\end{split}\]
元の問題に加えて、slack変数であるv0, w0が導入されました。
slack変数に付加される添字は、元の問題の添字を解析して自動的に決定されます。
また、目的関数にも v0 + w0 が加わっており、slack変数の値が小さくなるように最適化されます。
この問題で先ほどと同じように求解してみましょう。
目的関数に付加された bigM を別途インスタンスデータに定義する必要があります。
import numpy as np
from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter
instance_data = {
'A': np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
'b': np.array([10, 1, 20]),
'c': np.array([7, 8, 9]),
'bigM': 1000
}
interpreter = jm.Interpreter(instance_data)
ommx_instance = interpreter.eval_problem(replaced_problem)
# OMMXPySCIPOptAdapterを使用して問題を解く
solution = OMMXPySCIPOptAdapter.solve(ommx_instance)
solution.decision_variables
| kind | lower | upper | name | subscripts | description | substituted_value | value | |
|---|---|---|---|---|---|---|---|---|
| id | ||||||||
| 6 | integer | 0.0 | 100000000.0 | w0 | [0] | <NA> | <NA> | 0.0 |
| 7 | integer | 0.0 | 100000000.0 | w0 | [1] | <NA> | <NA> | 0.0 |
| 8 | integer | 0.0 | 100000000.0 | w0 | [2] | <NA> | <NA> | 0.0 |
| 0 | binary | 0.0 | 1.0 | x | [0] | <NA> | <NA> | 1.0 |
| 1 | binary | 0.0 | 1.0 | x | [1] | <NA> | <NA> | 1.0 |
| 2 | binary | 0.0 | 1.0 | x | [2] | <NA> | <NA> | 1.0 |
| 3 | integer | 0.0 | 100000000.0 | v0 | [0] | <NA> | <NA> | 4.0 |
| 4 | integer | 0.0 | 100000000.0 | v0 | [1] | <NA> | <NA> | 0.0 |
| 5 | integer | 0.0 | 100000000.0 | v0 | [2] | <NA> | <NA> | 0.0 |
これを見てみると、slack変数のうち、v0 の0番目の要素が0より大きい値を持っていることがわかります。
slack変数の値をなるべく小さくしようとしているにも関わらず、 v0 の0番目の要素が0より大きい値を持っていることから、元の問題において linear_constraint の i=0に対応する制約が実行不可能の原因である見当をつけることができます。
試しに、b[0]の値を小さくして、0番目の制約を緩和してみましょう。
import numpy as np
from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter
instance_data = {
'A': np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
'b': np.array([1, 1, 20]),
'c': np.array([7, 8, 9])
}
interpreter = jm.Interpreter(instance_data)
ommx_instance = interpreter.eval_problem(problem)
# OMMXPySCIPOptAdapterを使用して問題を解く
solution = OMMXPySCIPOptAdapter.solve(ommx_instance)
solution.decision_variables
| kind | lower | upper | name | subscripts | description | substituted_value | value | |
|---|---|---|---|---|---|---|---|---|
| id | ||||||||
| 0 | binary | 0.0 | 1.0 | x | [0] | <NA> | <NA> | 1.0 |
| 1 | binary | 0.0 | 1.0 | x | [1] | <NA> | <NA> | 1.0 |
| 2 | binary | 0.0 | 1.0 | x | [2] | <NA> | <NA> | 1.0 |
今度は実行可能解を得ることができました。 ついでに、slack変数を付加した問題も同じように説いてみると、slack変数の値が0になっていることがわかります。
import numpy as np
from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter
instance_data = {
'A': np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),
'b': np.array([1, 1, 20]),
'c': np.array([7, 8, 9]),
'bigM': 1000
}
interpreter = jm.Interpreter(instance_data)
ommx_instance = interpreter.eval_problem(replaced_problem)
# OMMXPySCIPOptAdapterを使用して問題を解く
solution = OMMXPySCIPOptAdapter.solve(ommx_instance)
solution.decision_variables
| kind | lower | upper | name | subscripts | description | substituted_value | value | |
|---|---|---|---|---|---|---|---|---|
| id | ||||||||
| 0 | binary | 0.0 | 1.0 | x | [0] | <NA> | <NA> | 1.0 |
| 1 | binary | 0.0 | 1.0 | x | [1] | <NA> | <NA> | 1.0 |
| 2 | binary | 0.0 | 1.0 | x | [2] | <NA> | <NA> | 1.0 |
| 3 | integer | 0.0 | 100000000.0 | v0 | [0] | <NA> | <NA> | 0.0 |
| 4 | integer | 0.0 | 100000000.0 | v0 | [1] | <NA> | <NA> | 0.0 |
| 5 | integer | 0.0 | 100000000.0 | v0 | [2] | <NA> | <NA> | 0.0 |
| 6 | integer | 0.0 | 100000000.0 | w0 | [0] | <NA> | <NA> | 0.0 |
| 7 | integer | 0.0 | 100000000.0 | w0 | [1] | <NA> | <NA> | 0.0 |
| 8 | integer | 0.0 | 100000000.0 | w0 | [2] | <NA> | <NA> | 0.0 |
このような方法で、slack変数を導入することで、実行不可能な問題を常に実行可能な問題に変換し、さらにslack変数の値を見てどの制約が悪さをしているのかを特定することができます。