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

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

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

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

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

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

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

定式化 弱い定式化

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

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

定式化 最大距離を最小化

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

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

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

実験結果

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

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

Miller-Tucker-Zemlinの定式化

上界と下界の変化 （８０点，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

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

対称な問題（DFJ定式化）

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

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

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

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

多品目ロットサイズ決定問題 段取り費用と在庫費用のトレードオフを最適化する多期間生産計画 多品目で共通の資源を使う容量制約付き問題は，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ベースのメタ解法を作ってみたが，同時間通常の探索を行う「打ち切り分枝限定法」の方が良い！

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

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

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

切断問題の定式化 初期パターンの線形緩和問題 切断パターンを使用する回数 [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) 初期切断パターン（どのアイテムを何回切断するか）

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

Pythonによる記述（主問題） K = len(t) 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) t: 初期パターンのリスト

Pythonによる記述（列生成） knapsack.update() while 1: 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")

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

グラフ彩色問題 解の対称性 点数 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秒

Fixed-K approach 悪い枝の数を最小化 枝の両端が同じ色（悪い枝）だと1

Fixed-K Approach +Binary Search 彩色数 Kの上界 UB，下界 LB while UB – LB >1: K= [ (UB+LB)/2 ] [ ] : Gaussの記号 if Fixd-K定式化の最適値が 0 then UB = K else LB = K

Improved Formulation Fixed-K Approach +Binary Search

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

多目的最適化 非劣解の集合をどのように生成するか？ （Gurobiは分枝限定法で見つかった複数の解を出力） 線形和によるスカラー化 制約に変換してスライドさせる方法 理想点からの距離を最小化する方法

非劣解 とスカラー化法 第二目的関数 この領域の解は Aに優越される スカラー化の 目的関数の方向 A スカラー化が見逃す非劣解 第一目的関数

二目的巡回セールスマン問題に対する制約スライド法の結果

制約スライド法と理想点法の比較

Model Gurobi Objects Variable Constraint GRBError addVar Column addSOS addConstr Constraint SOS Callbacks LinExpr QuadExpr

Model Object = Model(“name of model”） Methods Attributes Params ObjVal ModelSense Runtime Status SolCount ... AddVar AddConstr optimize update getVars relax ...

Prams (Parameters to control the optimizer) TimeLimit： limit of the computational time MIPGap：upper bound of the relative error MIPFocus：MIP search strategy SolutionNumber：number of solution LPMethod：solution method of linear optimization OutputFlag：0=Off，1=On ... E.g., : ModelObject.Params.TimeLimit=10

Model SCOP Objects Variable Linear Quadratic Alldiff addVariable(s) addConstraint Alldiff 詳細については http://www.logopt.com/scop/

Model Object = Model(“name of model”） Methods Attributes Params variables constraints AddVariable(s) AddConstraint optimize

Prams (Parameters to control the optimizer) TimeLimit： limit of the computational time Target：target penalty value RandomSeed: random seed number OutputFlag： 0=Off，1=On

Model OptSeq Objects Attribute Mode Resource Temporal addActivity addMode Model Mode addResource Resource addTemporal Temporal 詳細については http://www.logopt.com/OptSeq/

Model Object = Model(“name of model”） Methods Attributes Params activities modes resources temporals AddActivity AddResource AddTemporal optimize write

Prams (Parameters to control the optimizer) TimeLimit： limit of the computational time Makespan: True= makespan minimization, False=weighted tardiness minimization RandomSeed: random seed number OutputFlag： 0=Off，1=On