Python を用いた Gurobi 入門

1. Pythonを用いたGurobi入門 久保 幹雄

2. Why Python? • モジュールを読み込めば何でもできる！ • 最適化もできる！import gurobipy (MIP)import SCOP (CP) • グラフも描ける！import networkXimport matplotlib • 空を飛ぶことも！？import antigravity ? http://www.logopt.com/mikiokubo/に20分の入門ビデオ http://xkcd.com/353/

3. Python 1ページ解説 複合型 • リスト：任意の要素から成る順序型A=[] これだけで，スタック，キュー，連結リスト，ソート A.append(5), a=A.pop(), A.remove(“6”), A.sort() • 辞書： キーと 値の組から構成されるマップ型D= { }Hash関数と同じ．何でも辞書で高速に保管できる！ D[(1,2)] =6, D[“Hello”]=“こんにちは” 反復 for i in [1,2,5,6]: for key in D: print i*2 print key,D[key] => 2,4,10,12 =>キーと値を出力

4. k-median問題 • min-sum型の施設配置問題 • 顧客数 n=200，施設数 k=20（顧客上から選択） • Euclid距離，x,y座標は一様ランダム n=200, k=5 の例

5. 定式化 弱い定式化

6. Pythonでの実装（1） from gurobipy import * #gurobipyモジュールの読み込み # k-medianソルバーの関数 def solve(n,k,cost): model=Model(“median”) #モデルオブジェクトの生成 y={} #変数を表す辞書の準備 x={} 値 “127cm” 変数オブジェクト キー “Hanako”, (1,2) 写像

7. Pythonでの実装（２） for j in range(n): y[j]=model.addVar(obj=0,vtype=“B”,name="y"+str(j)) for i in range(n):x[i,j]=model.addVar(obj=cost[i,j], vtype=“B”,name="x"+str(i)+str(j)) model.update()

8. Pythonでの実装（３） for i in range(n): L=LinExpr() #線形表現オブジェクト L を空で作成 for j in range(n): L.addTerms(1,x[i,j]) #線形表現オブジェクトに項を追加 model.addConstr(lhs=L,sense= " = ",rhs=1,name="Demand"+str(i)) 線形表現（Linear Express） クラス LinExpr

9. Pythonでの実装（４） 中略 model.optimize() #最適化実行 print “Opt.value=”,model.ObjVal #目的関数値 edge=[] #選ばれた枝(i,j)のリスト for (i,j) in x: #変数xの辞書を反復 if x[i,j].X==1: #x[i,j]の解が1なら edge.append((i,j)) #リストedgeに(i,j)を追加 return edge

10. Pythonでの実装（５） import networkx as NX #networkXモジュールの読み込み import matplotlib.pyplot as P #描画の準備 P.ion() G = NX.Graph() #グラフオブジェクトGの生成 G.add_nodes_from(range(n)) #点の追加 for (i,j) in edge: #グラフに枝を追加 G.add_edge(i,j) NX.draw(G) #描画

11. 弱い定式化での結果n=200,k=20 Optimize a model with 401 Rows, 40200 Columns and 80400 NonZeros 中略 Explored 1445 nodes (63581 simplex iterations) in 67.08 seconds Thread count was 2 (of 2 available processors) Optimal solution found (tolerance 1.00e-04) Best objective 1.0180195861e+01, best bound 1.0179189780e+01, gap 0.0099% Opt.value= 10.1801958607

12. 上界と下界の変化

13. 強い定式化での結果 Optimize a model with 40401 Rows, 40200 Columns and 160400 NonZeros 中略 Explored 0 nodes (1697 simplex iterations) in 3.33 seconds(分枝しないで終了！） Thread count was 2 (of 2 available processors) Optimal solution found (tolerance 1.00e-04) Best objective 1.0180195861e+01, best bound 1.0180195861e+01, gap 0.0% Opt.value= 10.1801958607

14. 知見 • Big Mを用いない強い定式化が望ましい． • この程度の式なら，必要な式のみを切除平面として追加するような小細工は必要なし．（ただし，式の数が増え退化するので，大規模なＬＰを高速に解けるソルバーが前提）

15. 巡回セールスマン問題（ＴＳＰ） • すべての点をちょうど１回通る最短巡回路 • 切除平面法で85,900点の実際問題（対称ＴＳＰ）の最適解

16. Miller-Tucker-Zemlinの定式化

17. 上界と下界の変化（８０点，Euclid TSP） 強化した式 でないと．．． １日まわして Out of Memory!

18. 結果 Optimize a model with 6480 Rows, 6400 Columns and 37762 NonZeros 中略 Cutting planes: Gomory: 62 Implied bound: 470 MIR: 299 Zero half: 34 Explored 125799 nodes (2799697 simplex iterations) in 359.01 seconds Optimal solution found (tolerance 1.00e-04) Best objective 7.4532855108e+00, best bound 7.4525704995e+00, gap 0.0096% Opt.value= 7.45328551084

19. 知見 • 式を持ち上げ操作などで強化すると，高速化され，大きな問題例が解けるようになる． 　（そのためには多面体論の知識が多少必要なので，専門家に相談する）

20. 多品目ロットサイズ決定問題 • 段取り費用と在庫費用のトレードオフを最適化する多期間生産計画 • 多品目で共通の資源を使う容量制約付き問題は，MIPソルバーには難問（と言われてきた） • T=30期，P=24品目：Trigeiro, Thomas, McClain （1989年）の最大のベンチマーク

21. ロットサイズ決定問題標準定式化のフローモデルロットサイズ決定問題標準定式化のフローモデル 生産量(t) 在庫量(t) 在庫量(t-1) 期 t 需要量(t) 弱い定式化の原因 生産量(t)≦大きな数 “Large M” ×段取りの有無(t) 在庫量（t-1)+生産量(t)=需要量(t)+在庫量（t) 0-1変数

22. ロットサイズ決定問題施設配置定式化のフローモデルロットサイズ決定問題施設配置定式化のフローモデル s期に生産してt期まで在庫される量 期 s 期 t 需要量(t) s期に生産してt期まで在庫される量 ≦需要量(t)×段取りの有無(t) s期に生産してt期まで在庫される量 = 需要量(t)

23. 定式化のサイズと強さの比較 標準定式化 施設配置定式化 変数の数 変数の数 弱い定式化 強い定式化 =線形計画緩和が 整数多面体と一致 制約の数 制約の数 (S,l)不等式 切除平面 追加した 制約の数 n: 期数 強い定式化

24. 上界と下界の変化（標準定式化） 1800秒で最適解；これ以上大きな問題例は無理！

25. 上界と下界の変化（施設配置定式化） 40秒で最適解；T=100でも大丈夫！

26. 知見 • 従来では難問と言われてきたロットサイズ決定問題でもある程度までは大丈夫 • 緩和固定法（Federgruen, Meissner,Tzur “Progressive Interval Heuristics for Multi-Item Capacitated Lot-Sizing Problems” 2007）というMIPベースのメタ解法を作ってみたが，同時間通常の探索を行う「打ち切り分枝限定法」の方が良い！

27. k-center問題 • min-max型の施設配置問題 • 100顧客，10施設（顧客上から選択） • Euclid距離，x,y座標は一様ランダム

28. 定式化

29. 上界と下界の変化

30. 知見 • Min-max型の目的関数はMIPソルバーでは解きにくい（双対ギャップが大きい；上界も下界も悪いので，途中で止めても悪い解！）． • 例： Job shopスケジューリングの最大完了時刻（メイクスパン）最小化 • 制約計画としてモデル化して，上限を制約として扱うとうまくいく場合もある．

31. グラフ彩色問題 • 解の対称性 • 点数 n=40，彩色数上限 Kmax=10 • ランダムグラフ G(n,p=0.5) • 彩色数 8 が最適値

32. 定式化 弱い定式化

33. 上界と下界の変化（原定式化）点数 n=40，彩色数上限 Kmax=10 Optimize a model with 3820 Rows, 410 Columns and 11740 NonZeros Explored 17149 nodes (3425130 simplex iterations) in 1321.63 seconds

34. 定式化の改良 対称性の除去（番号の小さい方の色を優先して使う！） 特殊順序集合（SOS: Special Ordered Set） Type 1 （いずれか １つの変数が正になることの宣言）の追加 model.addSOS(1,変数リスト)

35. 上界と下界の変化（対称性除去） Optimize a model with 3829 Rows, 410 Columns and 11758 NonZeros Explored 4399 nodes (1013290 simplex iterations) in 384.53 seconds MIPFocus=2（最適性保証優先） で67秒，MIPFocus=3（下界優先） で70秒

36. 上界と下界の変化（+SOS） Optimize a model with 3829 Rows, 410 Columns and 11758 NonZerosExplored 109 nodes (58792 simplex iterations) in 22.02 seconds MIPFocus=2（最適性保証優先） で65秒，MIPFocus=3（下界優先） で126秒

37. 知見 • 解の対称性のある問題は，下界の改善がしにくいので，分枝限定法では解きにくい（モダンなソルバーは群論などを用いた工夫を入れてはいるが，あまり機能しない問題もある） • 対称性を除く工夫を入れると多少は改善 • SOS（Special Ordered Set）制約の宣言は損はない（ただし，探索手法の変更との相性は悪い．）

38. Gurobiクラス • モデルクラス Model モデルオブジェクト=Model(“モデル名”） • 変数クラス Var変数オブジェクト = モデルオブジェクト.addVar() • 制約クラス Constr制約オブジェクト = モデルオブジェクト.addConstr() • 線形制約クラス LinExpr線形制約オブジェクト=LinExpr() • 他にも，特殊順序集合クラス SOS，定数クラス GRB，コールバッククラス Callbacks，エラークラス GurobiError，列クラス Column がある．

39. 変数追加メソッド addVarの引数 • lb (下限) ：規定値 =0 • ub (上限) ：規定値=GRB.INFINITY（無限大） • obj (目的関数の係数)：＝０ • vtype (変数の種類): GRB.CONTINUOUS（連続変数； “C”でも可）, GRB.BINARY（2値変数；“B”）, GRB.INTEGER（整数変数；“I”）, GRB.SEMICONT（0もしくはある範囲内の連続変数；“S”でも可）, GRB.SEMIINT（0もしくはある範囲内の整数変数；“N”でも可） • name (名前): =“” • column (列): =なし（None型）

40. 制約追加メソッドaddConstrの引数 • ｌhs（左辺）：線形制約オブジェクトか定数 • sense（制約の向き）：GRB.LESS_EQUAL（≦；“<”でも可）, GRB.EQUAL（＝； “=”でも可）, GRB.GREATER_EQUAL （≧；“>”でも可） • rhs（右辺）：線形制約オブジェクトか定数 • name（名前）

41. モデルのその他の主要メソッド • optimize(コールバック関数）：最適化実行 • update：モデルの変更をGurobiに知らせる • getVars：変数オブジェクトを入れたリストを得る • relax：緩和問題を生成 • computeIIS：既約不整合部分系（実行不能になっている原因の制約と変数；Irreducible Inconsistent Subsystem：IIS）を計算 • addSOS（1 or 2, 変数リスト，重み（option））：特殊順序集合（Special Ordered Set：SOS）を追加

42. パラメータの設定 書式： setParam(“パラメータ名”,値） model.Params.パラメータ名 =値 例：混合整数計画（MIP）の相対誤差（終了条件）；規定値は 10-4 setParam(“MIPGap”,0.0) model.Params.MIPGap=0.0

43. よく使うパラメータ • TimeLimit：計算時間上限（秒） • MIPGap：相対誤差上限（規定値は 10-4） • MIPFocus：MIP探索戦略（0=バランス,1=暫定解改良,2=最適性の保証,3=限界値改良） • SolutionNumber：探索中に発見された解の番号 • LPMethod：線形計画の解法（0=主単体，1=双対単体，2=内点） • OutputFlag：出力制御（0=Off，1=On）

44. 属性（attributes） 書式： オブジェクト.setAttr(“属性名”,値） オブジェクト.属性名 =値 例：変数varの目的関数を100に変更 var.setAttr(“Obj”,100) var.Obj=100

45. よく使う属性（1） • モデル • ObjVal：最適目的関数値（最適化後のみ） • ModelSense：目的関数の方向（1=最小化，2=最大化） • Runtime：計算時間（最後の求解時の） • Status：モデルの状態（１から1３まで；2=OPTIMAL, 3=INFEASIBLE, 4=UNBOUNDED, 9=TIME_LIMIT, 11=INTERRUPTED, 13=SUBOPTIMAL) • SolCount：探索で見つかった解の数

46. よく使う属性（２） • 変数 • LB，UB：下限，上限 • Obj：目的関数の係数 • Vtype：変数の種類（“C”=連続，“B”=2値，“I”=整数，“S”＝半連続，“N”=半整数 • VarName：変数名 • X：解の値 • Xn：パラメータSolutionNumberで指定された番号の解の値 • RC：被約費用 • 制約 • Sense：制約の向き（“<”，“>”，“=”） • RHS：右辺定数 • ConstrName：制約名 • Pi：双対変数（連続最適化のみ） • Slack：余裕変数の値