Benders分解

Benders分解#

🌐 日本語 | English

概要#

Benders分解（ベンダース分解）は、バイナリ変数と連続変数が混ざった最適化問題（混合整数計画問題）を、2つの小さな問題に分けて交互に解く手法です。

マスター問題: バイナリ変数だけを決める問題
サブ問題: マスター問題が決めたバイナリ変数を固定したうえで、連続変数だけを決める問題

サブ問題を解くたびに、その結果からカット（マスター問題に追加する制約）を作ってマスター問題に戻します。カットが増えるほどマスター問題は元の問題に近づき、下界（最適値がこれより小さくならない値）と上界（実際に見つかった解の目的関数値）が両側から狭まっていきます。両者が一致したら、その解が最適解です。

少数のバイナリ変数の決定に多くの連続変数がぶら下がっている問題（例: どの拠点を開設するかを決め、開設した拠点での生産量を決める）で特に有効です。

JijZept Tools の BendersDecomposition クラスは、この一連の流れ（問題の分割・カットの生成・反復）をすべて自動で行います。ユーザーが用意するのは次の2つだけです。

JijModeling で記述した問題
マスター問題を解くソルバ（OMMX アダプタ形式の関数）

なお、連続変数のサブ問題は線形計画問題になるため、クラス内部で HiGHS（オープンソースの線形計画ソルバ）により自動的に解かれます。

例題#

小さな生産計画の問題を考えます。2つの拠点 A, B について「開設するかどうか」（バイナリ変数 \(x, y\)）と「開設した拠点での生産量」（連続変数 \(z_A, z_B\)）を決めます。開設には固定費（20, 30）、生産には単価（2, 3）がかかり、合計 5 の需要を満たす必要があります。

\[\begin{split} \begin{aligned} \min \quad & 20x + 30y + 2 z_A + 3 z_B \\ \mathrm{s.t.} \quad & x + y \le 1 \\ & z_A \le 10x, \quad z_B \le 10y \\ & z_A + z_B \ge 5 \\ & x, y \in \{0,1\}, \quad 0 \le z_A, z_B \le 10 \end{aligned} \end{split}\]

BendersDecomposition は、問題中の制約を変数の種類に応じて自動で分類します。

制約	含まれる変数	分類
`choose_one` (\(x + y \le 1\))	バイナリのみ	マスター問題
`capacity_A`, `capacity_B` (\(z_A \le 10x\) など)	バイナリ + 連続	結合制約（マスター問題とサブ問題をつなぐ制約）
`demand` (\(z_A + z_B \ge 5\))	連続のみ	サブ問題

import jijmodeling as jm

# 決定変数
x = jm.BinaryVar("x")  # 拠点 A を開設するか
y = jm.BinaryVar("y")  # 拠点 B を開設するか
z_a = jm.ContinuousVar("z_a", lower_bound=0, upper_bound=10)  # 拠点 A の生産量
z_b = jm.ContinuousVar("z_b", lower_bound=0, upper_bound=10)  # 拠点 B の生産量

problem = jm.Problem("facility_production")

# 目的関数: 開設固定費 + 生産コスト を最小化
problem += 20 * x + 30 * y + 2 * z_a + 3 * z_b

# 制約
problem += jm.Constraint("choose_one", x + y <= 1)     # 開設は1拠点まで（バイナリのみ）
problem += jm.Constraint("capacity_A", z_a <= 10 * x)  # 開設した拠点でのみ生産できる（結合）
problem += jm.Constraint("capacity_B", z_b <= 10 * y)  # 同上（結合）
problem += jm.Constraint("demand", z_a + z_b >= 5)     # 需要を満たす（連続のみ）

problem

\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{facility\_production} & & \\& & \min \quad \displaystyle 20 \cdot x + 30 \cdot y + 2 \cdot z_a + 3 \cdot z_b & \\\text{{s.t.}} & & & \\ & \text{capacity\_A} & \displaystyle z_a \leq 10 \cdot x & \\ & \text{capacity\_B} & \displaystyle z_b \leq 10 \cdot y & \\ & \text{choose\_one} & \displaystyle x + y \leq 1 & \\ & \text{demand} & \displaystyle z_a + z_b \geq 5 & \\\text{{where}} & & & \\& x & 0\text{-dim binary variable}\\& y & 0\text{-dim binary variable}\\& z_a & 0\text{-dim continuous variable}\\ & & \text{lower bound: }0 & \\ & & \text{upper bound: }10 & \\& z_b & 0\text{-dim continuous variable}\\ & & \text{lower bound: }0 & \\ & & \text{upper bound: }10 & \\\end{array}\end{split}\]

実行#

BendersDecomposition に問題とインスタンスデータ（jm.Placeholder に与える値の辞書。この例では Placeholder を使っていないため空）を渡してインスタンス化し、solve() を呼ぶだけです。

solve() の第1引数はマスター問題を解くソルバで、「OMMX の Instance を受け取り Solution（または SampleSet）を返す関数」であれば何でも構いません。ここでは HiGHS の OMMX アダプタをそのまま渡します。

verbose=True にすると、反復ごとに下界（LB）と上界（UB）が狭まっていく様子が表示されます。

from ommx_highs_adapter import OMMXHighsAdapter

from jijzepttools.modeling.algorithm.benders_decomposition import BendersDecomposition

benders = BendersDecomposition(problem, instance_data={})

result = benders.solve(
    OMMXHighsAdapter.solve,  # マスター問題を解くソルバ
    max_iterations=50,
    tolerance=1e-6,
    verbose=True,
)

Starting Benders decomposition with 2 binary vars, 2 continuous vars
Coupling constraints: 2
------------------------------------------------------------

Iteration   1: LB=-100000.000, UB=5010.000, Gap=105010.000, Feasible=True

Iteration   2: LB=-4970.000, UB=  30.000, Gap=5000.000, Feasible=True

Iteration   3: LB=  30.000, UB=  30.000, Gap=  -0.000, Feasible=True

Converged after 3 iterations!
Optimal value: 30.000000

結果の見方#

solve() は BendersResult を返します。主なフィールドは次のとおりです。

converged: 収束したか（下界と上界の差が tolerance 以下になったか）
iteration: 反復回数
lower_bound / upper_bound: 最終的な下界・上界
is_feasible: サブ問題が実行可能だったか
solution: 最終解（OMMX の Solution オブジェクト）

期待どおり、「固定費の安い拠点 A のみを開設し（\(x=1\)）、そこで需要分 5 を生産する」（目的関数値 \(20 + 2 \times 5 = 30\)）が得られます。

print("収束したか:", result.converged)
print("反復回数:", result.iteration)
print("目的関数値（上界）:", result.upper_bound)

# 最終解の変数値（OMMX Solution から取得）
result.solution.decision_variables_df[["name", "value"]]

収束したか: True
反復回数: 3
目的関数値（上界）: 30.0

	name	value
id
0	x	1.0
1	y	0.0
2	z_a	5.0
3	z_b	0.0

内部の動きを覗く#

インスタンス化した時点で、変数と制約は自動的にマスター問題側・サブ問題側・結合制約に分類されています。また、反復のたびに生成されたカットは cuts に蓄積されます。

1反復ずつ手動で進めたい場合は、solve() の代わりに iterate()（1反復だけ実行して BendersResult を返す）を繰り返し呼ぶこともできます。

print("マスター問題側のバイナリ変数:", [v.name for v in benders.binary_vars])
print("サブ問題側の連続変数:", [v.name for v in benders.continuous_vars])
print("結合制約:", [c.name for c in benders.coupling_constraints])
print("追加されたカットの数:", len(benders.cuts))

マスター問題側のバイナリ変数: ['x', 'y']
サブ問題側の連続変数: ['z_a', 'z_b']
結合制約: ['capacity_A', 'capacity_B']
追加されたカットの数: 2

マスター問題のソルバを差し替える#

マスター問題のソルバは「OMMX の Instance を受け取って解く関数」なら何でもよいので、他のソルバへの差し替えは1行で済みます。たとえば SCIP（オープンソースの混合整数計画ソルバ）の OMMX アダプタを使う場合は次のとおりです。

from ommx_pyscipopt_adapter import OMMXPySCIPOptAdapter

benders_scip = BendersDecomposition(problem, instance_data={})
result_scip = benders_scip.solve(OMMXPySCIPOptAdapter.solve, max_iterations=50, tolerance=1e-6)

print("収束したか:", result_scip.converged)
print("目的関数値:", result_scip.upper_bound)

収束したか: True
目的関数値: 30.0

マスター問題はバイナリ変数だけの問題なので、アニーリング等の QUBO 系サンプラー（例: OpenJij）をマスター問題のソルバとして使うこともできます。サンプラーが SampleSet（複数の解の集まり）を返す場合はマルチカットが有効になり、1回の反復で複数のサンプルからまとめてカットを生成できるため、反復回数を減らせることがあります。

主なオプション#

BendersDecomposition のコンストラクタでは、次のようなオプションを指定できます（詳細は各引数の docstring を参照）。

slack_M: サブ問題が実行不能なときでも計算を続けられるように導入されるスラック変数（制約からの逸脱を表す補助変数）へのペナルティ係数（既定値 1000）
use_two_phase: 実行可能性の確保（フェーズ1）と最適化（フェーズ2）を分けて行う2段階方式を使うか
normalize_cuts: カットの係数を正規化するか（係数の桁が大きく異なる場合に有効）
max_cuts: 保持するカット数の上限（超過時は活性でないカットから削除。既定値 100）
max_samples: SampleSet を受け取ったとき、カット生成に使うサンプル数の上限（目的関数値のよい順に選択）