Benders Decomposition(ベンダー分解)は問題を二つのmain, sub問題に分割することで問題を解きやすくするテクニックです。次のような性質がある問題に有効です。Kohji Nishimuraさんの「双対変数を直感的に理解したい」から引用しています。
決定変数にyが一部含まれているが、もしこのyがなかったら解くのが簡単な問題なのになぁ、といったシチュエーションは数理最適化にて頻繁に起こる。例えば
- yがなければ残った問題が複数の小さな問題に分解できる
- yが離散変数 (バイナリ変数等)であるため、yさえ固定できてしまえば残りの問題が線形計画問題として解ける
双対変数はややこしいイメージがありましたが、上記の記事で分かりやすく説明されており、理解が深まりました。 今回は簡単な問題に対して、Benders DecompositionをGurobipyを用いて実装したという記事です。 Benders Decompositionについては次のsnowberryfieldさんの記事「Benders分解法」が詳しいです。
テスト問題
下記のサイトの問題を用いました。そのまま解くと最適な目的関数値は2。
https://homepages.rpi.edu/~mitchj/matp6620/benders_eg/benders_eg_slides.html
from gurobipy import Model, GRB # Create a new model m = Model("mip_problem") # Create variables x = m.addVars(4, vtype=GRB.CONTINUOUS, name="x") y = m.addVars(2, vtype=GRB.BINARY, name="y") # Set objective m.setObjective(8*x[0] + 9*x[1] + 5*x[2] + 6*x[3] - 15*y[0] - 10*y[1], GRB.MAXIMIZE) # Add constraints m.addConstr(x[0] + x[2] <= 1) m.addConstr(x[0] + x[3] <= 1) m.addConstr(x[1] + x[2] <= 1) m.addConstr(x[1] + x[3] <= 1) m.addConstr(x[0] - y[0] <= 0) m.addConstr(x[1] - y[0] <= 0) m.addConstr(x[2] - y[1] <= 0) m.addConstr(x[3] - y[1] <= 0) m.addConstr(-y[0] - y[1] <= -1)
gurobipy実装
上記のブログに記載の手順とは異なりますが、行なっている処理は同じです（本稿の方が多少無駄が含まれます）。 テスト問題のバイナリ変数yがsubから分離する決定変数たちです。 上記のブログではsubで双対問題を解いていますが、下のコードでは対応する主問題を解いてその後双対変数の値を取得することで同様の制約をmain問題に追加しています。
main問題の初期化
import math from gurobipy import Model, GRB # Creat main problem main = Model("main") main_y = main.addVars(2, vtype=GRB.BINARY, name="y") main_t = main.addVar(lb=-GRB.INFINITY, vtype=GRB.CONTINUOUS, name="t") main.addConstr(- main_y[0] - main_y[1] <= -1) main.setObjective(main_t) main.Params.OutputFlag = 0 # silent/verbose mode
sub問題の初期化
# Create sub problem sub = Model("sub") x = sub.addVars(4, vtype=GRB.CONTINUOUS, name="x") y = sub.addVars(2, lb=0, ub=1, vtype=GRB.CONTINUOUS, name="y") sub.addConstr(x[0] + x[2] <= 1) sub.addConstr(x[0] + x[3] <= 1) sub.addConstr(x[1] + x[2] <= 1) sub.addConstr(x[1] + x[3] <= 1) sub.addConstr(x[0] - y[0] <= 0) sub.addConstr(x[1] - y[0] <= 0) sub.addConstr(x[2] - y[1] <= 0) sub.addConstr(x[3] - y[1] <= 0) d0 = sub.addConstr(0 <= 0) # dummy d1 = sub.addConstr(0 <= 0) # dummy sub.setObjective(-8*x[0] - 9*x[1] - 5*x[2] - 6*x[3] + 15*y[0] + 10*y[1]) sub.Params.OutputFlag = 0 # silent/verbose mode
Benders分解ループ
# Set initial solution randomly opt_y = [0 for i in range(2)] main_obj = float("-inf") # Benders decomposition index = 0 while True: # Update sub problem sub.remove([d0, d1]) d0 = sub.addConstr(y[0] == opt_y[0]) d1 = sub.addConstr(y[1] == opt_y[1]) sub.update() # Solve sub problem sub.optimize() sub_obj = sub.ObjVal l = [d0.Pi, d1.Pi] # get dual variable values # Output logging print(f"Round {index}") print(f" opt_y = {opt_y}") print(f" main_obj = {main_obj}") print(f" sub_obj = {sub_obj}") print(f" lambda = {l}") # Check termination condition if math.isclose(main_obj, sub_obj, abs_tol=1e-6): break # Add constraint to the main problem c = main.addConstr(main_t >= sub_obj + quicksum(l[i] * (main_y[i] - opt_y[i]) for i in range(2))) main.update() print(f" add {main_t >= sub_obj + quicksum(l[i] * (main_y[i] - opt_y[i]) for i in range(2))}") # Solve main problem main.optimize() main_obj = main.ObjVal opt_y = [main_y[i].X for i in range(2)] index += 1 print() print("Result") for v in sub.getVars(): print(f' {v.varName} = {v.x}')
以下が求解時のログ。最適化の方向を最小化にしているため符号が逆転していますが、オリジナルの問題を直接解いた場合の目的関数値2と一致していますね。
Round 0
opt_y = {0: 0, 1: 0}
main_obj = -inf
sub_obj = 0.0
lambda = [-2.0, -1.0]
add <gurobi.TempConstr: t >= -2.0 y[0] + -1.0 y[1]>
Round 1
opt_y = {0: 1.0, 1: 1.0}
main_obj = -3.0
sub_obj = 8.0
lambda = [15.0, 10.0]
add <gurobi.TempConstr: t >= -17.0 + 15.0 y[0] + 10.0 y[1]>
Round 2
opt_y = {0: 1.0, 1: 0.0}
main_obj = -2.0
sub_obj = -2.0
lambda = [0.0, -1.0]
add <gurobi.TempConstr: t >= -2.0 + 0.0 y[0] + -1.0 y[1]>
Result
x[0] = 1.0
x[1] = 1.0
x[2] = 0.0
x[3] = 0.0
y[0] = 1.0
y[1] = 0.0
