サプライ・チェインにおける様々な最適化問題を解くための 統一言語 東京海洋大学 久保幹雄

目次 Pythonを使えば最適化も簡単！ サプライ・チェイン最適化のお話 Python入門 色々な問題を解いてみた＋多少の知見 統一モデル

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

Python 1ページ解説 複合型 反復 for i in [1,2,5,6]: for key in D: リスト：任意の要素から成る順序型 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 =>キーと値を出力

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

定式化 弱い定式化

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

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

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

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

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

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

上界と下界の変化

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

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

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

Miller-Tucker-Zemlinの定式化

AMPLでの実装（１） param n >=0; set V := 1..n ; #点集合 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];

AMPLでの実装（２） Degree1 {i in V}: sum {(i,j) in A } x[i,j] =1 ; #出次数制約 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]; #持ち上げ上界制約

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

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

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

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

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

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

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

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

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

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

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

定式化

上界と下界の変化

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

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

定式化 弱い定式化

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)) y[k]=model.addVar(obj=1,vtype=“B”,name="y"+str(k)) model.update()

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

上界と下界の変化（原定式化） 点数 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

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

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

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

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