Python を用いた お気楽最適化とその実践

1. Pythonを用いたお気楽最適化とその実践 久保 幹雄

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

3. Pythonをお薦めする訳 キーワード（覚えるべき予約後）が30程度と少ない． 字下げの強要で，誰でも読みやすいプログラム 短時間で開発可能（行数が短く，モジュール豊富） 変数の宣言必要なし インタープリタ（コンパイルする必要なし；並列用コンパイラもついている） メモリ管理も必要なし 多くのプラットフォームで動作（Windows, Mac, Linux) オブジェクト指向（すべてがオブジェクト） しかもフリーソフト

4. 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 =>キーと値を出力

5. How to solve hard optimization problems quickly by Python Metaheuristics Python Source Codes Developed by: Me & J. Pedro Pedroso MIP solver Gurobi Developed by: Zonghao Gu, Edward Rothberg，Robert Bixby VS. ? + Free academic license + Python Interface

6. MIP, CP, Scheduling Solvers • 混合整数計画(MIP)ソルバー: Gurobi • 制約計画(CP)ソルバー: SCOP (Solver for COnstraint Programming) • スケジューリング問題に特化したソルバー: OptSeq II モデリング(定式化）のコツ使い分けのコツ 「新時代の最適化-Python言語を用いた最適化入門-」近代科学社（執筆中）参照

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

8. 定式化 弱い定式化

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

10. 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()

11. 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

12. 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

13. 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) #描画

14. 弱い定式化での結果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

15. 上界と下界の変化

16. 強い定式化での結果 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

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

18. k-center問題 • min-max型の施設配置問題 • 施設は顧客上から選択 • Euclid距離，x,y座標は一様ランダム Which is the solution of the k-center problem?

19. 定式化 最大距離を最小化

20. 上界と下界の変化 100顧客，10施設　（k-medianは200施設，20顧客で3秒）

21. k-被覆問題としての定式化 被覆されない顧客数最小化 被覆されなければ１の変数 距離がθ以下のとき 1のパラメータ

22. k-被覆問題+Binary Search 顧客・施設間の距離の上限 UB，下限 LB while UB – LB >ε: θ= (UB+LB)/2 if K-被覆問題の最適値が 0 then UB = θ else LB = θ

23. 実験結果

24. 知見 • Min-max型の目的関数はMIPソルバーでは解きにくい（双対ギャップが大きい；上界も下界も悪いので，途中で止めても悪い解！）． • 例： Job shopスケジューリングの最大完了時刻（メイクスパン）最小化 =>スケジューリング問題に特化したソルバーの方が良い！ • k-被覆問題を用いた二分探索だとある程度min-max型の目的関数を回避できる． • min-max型の問題は，制約計画としてモデル化して，上限を制約として扱うと良い．

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

26. Miller-Tucker-Zemlinの定式化

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

28. 結果 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

29. 単一フロー定式化 多少強化

30. 対称な問題（DFJ定式化）

31. 部分巡回路除去制約を追加線形緩和問題

32. 部分巡回路除去制約を追加整数解 連結成分を求めて切除平面を追加

33. 部分巡回路除去制約を追加最適解

34. 知見 • 式を持ち上げ操作などで強化すると，高速化され，大きな問題例が解けるようになる． 　（そのためには多面体論の知識が多少必要なので，専門家に相談．ちょっと式をいじるだけで劇的に改善する場合もある！） • 他にも単品種フロー，多品種フロー定式化などがある．（下界が強い方が良いとは限らない！LPのサイズと解きやすさ（退化の具合）を考慮して定式化を選択！） • 対称な問題なら分枝カット法が良い．

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

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

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

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

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

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

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

42. ビンパッキング問題 アイテムj=1,…,n のサイズsj サイズBのビン

43. 切断問題 データが離散値をとるビンパッキング問題と同じ

44. ビンパッキング問題の定式化 ビンjを使うとき１ アイテムiをビンjに入れるとき１ 強い定式化の ための追加式 解の対称性のためあまり大きい問題例は解けない！

45. 切断問題の定式化初期パターンの線形緩和問題切断問題の定式化初期パターンの線形緩和問題 切断パターンを使用する回数 [4, 0, 0, 0, 0, 0, 0] [0, 3, 0, 0, 0, 0, 0] [0, 0, 2, 0, 0, 0, 0] [0, 0, 0, 1, 0, 0, 0] [0, 0, 0, 0, 1, 0, 0] [0, 0, 0, 0, 0, 1, 0] [0, 0, 0, 0, 0, 0, 1] 最適双対変数 (1/4,1/3,1/2,1,1,1,1) 初期切断パターン（どのアイテムを何回切断するか）

46. 新しいパターンの生成ナップサック問題 y=(2, 0, 0, 1, 0, 0, 0)　　最適値 1.5 列を追加して 再び求解

47. Pythonによる記述（主問題） K = len(t) master = Model("MP") x = {} for k in range(K): x[k] = master.addVar(obj=1, vtype="I", name="x[%d]"%k) master.update() orders={} for i in range(m): coef = [t[k][i] for k in range(K) if t[k][i] > 0] var = [x[k] for k in range(K) if t[k][i] > 0] orders[i] = master.addConstr(LinExpr(coef,var), ">", q[i], name="Order[%d]"%i) master.update() t: 初期パターンのリスト

48. Pythonによる記述（列生成） knapsack.update() knapsack.optimize() if knapsack.ObjVal < 1+EPS: break col = Column() for i in range(m): if t[K][i] > 0: col.addTerms(t[K][i], orders[i]) x[K] = master.addVar(obj=1, vtype="I", name="x[%d]"%K, column=col) master.update() K += 1 master.optimize() while 1: iter += 1 relax = master.relax() relax.optimize() pi = [c.Pi for c in relax.getConstrs()] knapsack = Model("KP") knapsack.ModelSense=-1 y = {} for i in range(m): y[i] = knapsack.addVar(obj=pi[i], ub=q[i], vtype="I", name="y[%d]"%i) knapsack.update() L = LinExpr(w, [y[i] for i in range(m)]) knapsack.addConstr(L, "<", B, name="width")

49. 列生成法の応用事例（タンカースケジューリング問題）列生成法の応用事例（タンカースケジューリング問題） • タンカーのパーティションを考慮した複雑な実際問題 • 海技研からの委託 • 可能なパターンの列挙，集合被覆問題による選択を2回反復した解法 • 従来法（CP）：9時間かけて実行不能 • 開発法：10分程度で求解

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