Python を用いた Gurobi 入門 久保 幹雄. Why Python? モジュールを読み込めば何 でもできる！ 最適化もできる！ import gurobipy (MIP) import SCOP (CP) グラフも描ける！ import networkX import matplotlib.

Presentation on theme: "Python を用いた Gurobi 入門 久保 幹雄. Why Python? モジュールを読み込めば何 でもできる！ 最適化もできる！ import gurobipy (MIP) import SCOP (CP) グラフも描ける！ import networkX import matplotlib."— Presentation transcript:

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

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

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={} キー “Hanako”, (1,2) 写像 値 “127cm” 変数オブジェクト

7 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() Python での実装（２）

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 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 弱い定式化での結果 n=200,k=20

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 AMPL での実装（１） param n >=0; set V := 1..n ; # 点集合 set V0 := 2..n; # 出発地点 1 以外の点集合 set A :=V cross V; # 枝集合 = 点集合の直積（ 2 つ組） param c { A } >= 0; # 枝の距離 var x { A } binary ; # 枝を使うとき 1 ，それ以外のとき 0 の 0-1 変数 var u { V0 } >=1,<=n-1; # 点のポテンシャル minimize total_cost: sum {(i,j) in A} c[i,j] * x[i,j];

18 AMPL での実装（２） Degree1 {i in V}: sum {(i,j) in A } x[i,j] =1 ; # 出次数制約 Degree2 {i in V}: sum {(j,i) in A } x[j,i] =1 ; # 入次数制約 MTZ{ (i,j) in A: i != j and j!=1 and i!=1}: u[i]+1 -(n-1)*(1-x[i,j]) + (n-3)*x[j,i]<=u[j]; # 持ち上げ MTZ 制約 LiftedLB{ i in V0}: 1+(1-x[1,i]) +(n-3)*x[i,1] <= u[i]; # 持ち上げ下界制約 LiftedUB{ i in V0}: u[i] <=(n-1)-(1-x[i,1])-(n-3)*x[1,i]; # 持ち上げ上界制約

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

20 結果 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

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

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

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

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

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

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

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

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

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

30 定式化

31 上界と下界の変化

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

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

34 定式化 弱い定式化

35 Python での実装（１） from gurobipy import * model=Model("gcp") x={} y={} for i in range(n): for k in range(K): x[i,k]=model.addVar(obj=0, vtype="B",name="x"+str(i)+str(k)) for k in range(K): y[k]=model.addVar(obj=1,vtype=“B”,name="y"+str(k)) model.update()

36 Python での実装（２） for i in range(n): L=LinExpr() for k in range(K): L.addTerms(1,x[i,k]) model.addConstr(lhs=L,sense= " = ",rhs=1,name="const"+str(i)) 中略 model.optimize() print "Opt.value=",model.ObjVal for v in model.getVars(): if v.X>0.001: print v.VarName,v.X

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

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

39 上界と下界の変化（対称性除 去） 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 秒

40 上界と下界の変化（ +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 秒

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

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

43 変数追加メソッド 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 型）

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

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

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

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

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

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

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