この記事について
これまでの3記事で「LP緩和と上界・下界」「分枝限定法」「切除平面法」を説明してきました。今回はこれらを統合した 分枝カット法(Branch and Cut) と、実際の整数計画ソルバーが内部で何をしているのかを説明します。
「なんとなく Gurobi を使っているけど、計算が終わらないときどうすればいいの?」という疑問への答えにもなるかと思います。
※もちろん実際のソルバー開発はもっといろいろ考えて作られているかと思うので、あくまで一要素かなと思います。
分枝カット法(Branch and Cut)とは
分枝カット法は分枝限定法と切除平面法を組み合わせたアルゴリズムで、現代の整数計画ソルバーのほぼすべてが採用している手法です。
分枝限定法だけ: 問題を分割 → LP を解く → 枝刈り → ...
切除平面法だけ: LP にカットを追加 → LP を解く → カットを追加 → ...
分枝カット法: 上記を組み合わせ
基本的な流れ
① LP 緩和を解く
②-A LP 最適解が整数解 → 暫定値を更新(終端)
②-B LP が実行不能 → 枝刈り(終端)
②-C LP 緩和が暫定値を改善できない → 枝刈り(終端)
③ カット生成を試みる
→ カットが見つかった → LP にカットを追加して ① に戻る
→ カットが見つからない → ④へ
④ 分枝(変数固定して部分問題を生成)→ ① に戻る
各ノードで「まずカットを生成して LP 緩和を改善」してから「分枝が必要なら分枝」という戦略で、探索木のサイズを大幅に削減できます。
分枝カット法のイメージ
根ノード(LP緩和)
↓ カット追加 × 数回 → LP 緩和値が改善
↓ まだ整数解でない
↓ 分枝 → 部分問題 P₁, P₂ を生成
P₁: カット追加 → また分枝 → P₃, P₄ ...
P₂: カット追加 → 整数解が見つかる
実際のソルバーで使われるカット
Gurobi・CPLEX などの商用ソルバーは Gomory カット以外にも、問題構造を利用した強力なカットを複数種類使います。
| カット名 | 概要 |
|---|---|
| Gomory Fractional Cut | LP 最適辞書から自動生成(汎用的) |
| Mixed Integer Rounding (MIR) Cut | 混合整数制約から導出 |
| Clique Cut | 0-1 変数の相互排除関係を利用 |
| Cover Cut | ナップサック的な構造を利用 |
| Flow Cover Cut | ネットワーク構造を利用 |
| Lift-and-Project Cut | 複数変数の積の線形化 |
これらのカットは Gurobi の Cuts パラメータで制御できます。
Gurobi で分枝カット法の動きを観察する
Gurobi のログには分枝カット法の状態がリアルタイムで表示されます。
import gurobipy as gp
from gurobipy import GRB
import numpy as np
# ランダムな MIP を生成して解く
np.random.seed(42)
n = 30 # 変数数
m = 15 # 制約数
c = np.random.randint(1, 20, n).tolist()
A = np.random.randint(0, 10, (m, n)).tolist()
b = np.random.randint(30, 60, m).tolist()
model = gp.Model("random_mip")
model.setParam('OutputFlag', 1) # ログを表示
x = model.addVars(n, vtype=GRB.BINARY, name="x")
model.setObjective(gp.quicksum(c[j]*x[j] for j in range(n)), GRB.MAXIMIZE)
for i in range(m):
model.addConstr(gp.quicksum(A[i][j]*x[j] for j in range(n)) <= b[i])
model.optimize()
print(f"\n最適値: {model.ObjVal:.1f}")
print(f"MIP Gap: {model.MIPGap*100:.2f}%")
print(f"探索ノード数: {model.NodeCount:.0f}")
出力例(Gurobi ログの一部):
Gurobi 11.0.0 ...
...
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 237.00 0 15 - 237.00 - - 0s
0 0 235.12 0 10 - 235.12 - - 0s ← カット追加
0 0 228.10 0 8 210.00 228.10 8.62% - 0s ← 整数解発見
...
Optimal solution found (tolerance 1.00e-04)
Best objective 2.280000000000e+02, best bound 2.280000000000e+02, gap 0.0000%
ログの見方:
- Nodes: 探索済みノード数
- Incumbent: 現在の暫定解の目的関数値(これまでに見つかった最良の整数解の値)
- BestBd: 現在の LP 緩和の最適値(理論上の限界値)
- Gap: MIP ギャップ(Incumbent と BestBd の相対差)
最大化問題なら Incumbent が下界・BestBd が上界、最小化問題ならその逆になります。
計算が終わらないときの原因と対処
原因1:LP 緩和の求解自体が遅い
LP 緩和(整数条件を外した問題)を解くだけでも計算時間がかかる場合。
対処法
- 変数の数・制約の数を削減(不要な変数を事前に削除)
- 数値の桁数を揃える(係数の大小差が激しいと数値的に不安定)
原因2:LP 緩和の値の改善が遅い(LP 緩和が弱い)
LP 緩和の最適値が MIP 最適値から遠く離れており、カットが効かない場合。
対処法
- より強いカットを使う(Gurobi のカット設定を変更)
- 問題を再定式化して LP 緩和をタイトにする
- 対称性の除去(同等の解が多数ある場合)
原因3:暫定解の改善が遅い(良い整数解が見つからない)
分枝してもなかなか実行可能整数解が見つからない場合。
対処法
-
ヒューリスティクスで初期整数解を与える(
model.setStart()) - Gurobi の Heuristics パラメータでヒューリスティクス頻度を上げる
実践的なパラメータ設定
import gurobipy as gp
model = gp.Model()
# ... 問題定義 ...
# 時間制限(秒)
model.setParam('TimeLimit', 300)
# MIP ギャップ(1% 以内で終了)
model.setParam('MIPGap', 0.01)
# カット生成の積極性(-1=自動, 0=なし, 1=控えめ, 2=積極的, 3=非常に積極的)
model.setParam('Cuts', 2)
# ヒューリスティクスの時間配分(0〜1、デフォルト0.05)
model.setParam('Heuristics', 0.2)
# 深さ優先 vs 最良優先(0=自動, 1=深さ優先, 2=最良優先)
model.setParam('NodeMethod', 1)
model.optimize()
# 結果確認
if model.Status == gp.GRB.OPTIMAL:
print(f"最適解が見つかりました")
elif model.Status == gp.GRB.TIME_LIMIT:
print(f"時間制限で終了(Gap: {model.MIPGap*100:.2f}%)")
elif model.Status == gp.GRB.INFEASIBLE:
print(f"実行不能")
print(f"暫定解の値(Incumbent): {model.ObjVal:.2f}")
print(f"LP 緩和の値(BestBd): {model.ObjBound:.2f}")
print(f"MIP Gap: {model.MIPGap*100:.4f}%")
print(f"探索ノード数: {model.NodeCount:.0f}")
上界と下界の時間推移を可視化
ソルバーの収束過程を可視化すると、どこで詰まっているかが見えてきます。
import gurobipy as gp
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Yu Gothic'
upper_bounds = []
lower_bounds = []
times = []
def my_callback(model, where):
if where == gp.GRB.Callback.MIP:
obj = model.cbGet(gp.GRB.Callback.MIP_OBJBST) # 上界値
bnd = model.cbGet(gp.GRB.Callback.MIP_OBJBND) # 下界値
t = model.cbGet(gp.GRB.Callback.RUNTIME)
upper_bounds.append(obj)
lower_bounds.append(bnd)
times.append(t)
# model.optimize(my_callback) # コールバックを渡して最適化
# 収束曲線を描画
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(times, upper_bounds, label='上界値(暫定解)', color='#e74c3c', linewidth=1.5)
ax.plot(times, lower_bounds, label='下界値(LP緩和)', color='#3498db', linewidth=1.5)
ax.fill_between(times, lower_bounds, upper_bounds, alpha=0.15, color='#95a5a6')
ax.set_xlabel('計算時間(秒)')
ax.set_ylabel('目的関数値')
ax.set_title('分枝カット法の収束曲線')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('branch_cut_convergence.png', dpi=150)
上界値と下界値が収束していく様子が見える図が出力されます。
- 上界値(暫定解):ステップ状に減少する
- 下界値(LP緩和):滑らかに増加する
- ギャップ:時間とともに縮小し、一致した時点で最適解
分枝限定法・切除平面法・分枝カット法の関係
| 手法 | 下界の改善方法 | 特徴 |
|---|---|---|
| 分枝限定法 | 部分問題ごとに LP を解く | 汎用的。カットなし |
| 切除平面法 | LP にカットを追加 | カット次第で強力。収束が遅いことも |
| 分枝カット法 | 部分問題ごとにカット+LP | 両者のいいとこ取り。現代ソルバーの標準 |
Gurobi が速い理由:
- 強いカットを複数種類自動生成してLP緩和をタイト化
- ヒューリスティクスで良い整数解を早期に発見して上界を下げる
- 問題構造を認識したプリソルブ(前処理)で変数・制約を削減
- これらすべてを並列計算で高速化
まとめ
- 分枝カット法は分枝限定法と切除平面法を組み合わせた現代ソルバーの標準手法
- 各ノードで「まずカットを試みて下界を改善、それでも整数解が出なければ分枝」
- 計算が終わらないときは「LP緩和が弱い」「上界が得られていない」が主な原因
-
TimeLimitで時間制限、MIPGapで許容精度を設定することが実践上重要 - ソルバーのログや MIP Gap を確認することで、解の品質と収束状況を把握できる