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

Slides:



Advertisements
Similar presentations
組合せ最適化輪講 2.3 連結性 川原 純. 2.3 連結性 内容 – グラフ上の節点をすべてたどるアルゴリズム 計算機上でのグラフの表現 – 強連結成分を求めるアルゴリズム トポロジカル順序を求める方法も – k- 連結、 k- 辺連結について – 2- 連結グラフの耳分解について.
Advertisements

図示、可視化モジュール ~ pylab と numpy を ちょっと~. pylab とは? ・数学や統計的なグラフを生成するモ ジュール ・インストール pip や easy install からのインストールを推奨 →numpy モジュールなどの前提としている。 Anaconda の場合は標準.
情報基盤アルゴリズムとして のメタヒューリスティクスの 研究 茨木俊秀、巳波弘佳、藤原洋志、 千葉英史、関口良行(関西学院大 学)、 藤重悟(京都大学)、 柳浦睦憲(名古屋大学)、 野々部宏司(法政大学)、 梅谷俊治(電気通信大学)
模擬国内予選2013 Problem F テトラ姫のパズル 原案:須藤 解答:大友、須藤 解説:須藤.
2 dimensional bit vector approach Mikio Kubo. TIGER/Line graph DC.tmp 9559 nodes and arcs Shortest Paths between all border nodes of 2 regions.
Pythonを用いた お気楽最適化とその実践
制約付き最短路問題に対する実験的解析 上智大学 宮本裕一郎 o 1 d (7, 15) (9,4) (12,4)
NetworkXによるネットワーク分析 東京海洋大学 久保 幹雄.
JavaScript プログラミング入門 2006/11/10 神津.
ネットワーク理論講義補助資料 Text. 組合せ最適化とアルゴリズム 4.5 節 主・双対法 pp
「基礎OR」/「OR演習」 第2回 10/06/2009 森戸 晋.
ラベル付き区間グラフを列挙するBDDとその応用
プログラミング言語としてのR 情報知能学科 白井 英俊.
近似アルゴリズム 第10章 終了時刻最小化スケジューリング
サプライ・チェインにおける様々な最適化問題を解くための 統一言語
数理最適化の応用例と 実験的解析 東京海洋大学 久保 幹雄.
二千十三年五月 あたらしい数理最適化 出版記念セミナー 主催 近代科学社 オクトーバースカイ 共催 構造計画研究所
スレッドの同期と、スレッドの使用例 スレッドの同期 Lockオブジェクト: lockオブジェクトの生成
2017/3/10 スケジューリング最適化 東京海洋大学 久保 幹雄.
最短路問題をGurobiで解こう! 流通最適化工学 補助資料.
アルゴリズムとプログラミング (Algorithms and Programming)
Lightweight Language Weekend ls-lRシェル
整数計画法を用いた スリザーリンクの解法 杉村 由花 (東京大学)
最適化ソルバーのための Python言語入門
地理情報システム論 第3回 コンピュータシステムおける データ表現(1)
モード付き並列機械における オンラインスケジューリング
遺伝アルゴリズムによる NQueen解法 ~遺伝補修飾を用いた解探索の性能評価~
土木計画学 第11回(12月21日) 土木計画と説明責任 計画における代替案の作成1 担当:榊原 弘之.
Bottle/Pythonによる Webアプリ入門
「基礎OR」/「OR演習」 第3回 10/13/2009 森戸 晋.
整数計画法を用いた ペグソリティアの解法 ver. 2.1
データ構造とアルゴリズム 分割統治 ~ マージソート~.
1章前半.
Probabilistic Method 6-3,4
第7章 データベース管理システム 7.1 データベース管理システムの概要 7.2 データベースの格納方式 7.3 問合せ処理.
ML 演習 第 7 回 新井淳也、中村宇佑、前田俊行 2011/05/31.
ピカチュウによる オブジェクト指向入門 (新版)
第7回 条件による繰り返し.
情報工学総合演習 D-I 近似アルゴリズム 埼玉大学 理工学研究科 山田 敏規、 橋口 博樹、 堀山 貴史
割当て問題 • 割当問題の記法・定式化 • 拡張 • 特殊ケース(マッチング) • 3種類のものを割当てる問題.
大阪市立大学 学術情報総合センター 大西克実
ネットワーク理論 Text. Part 3 pp 最短路問題 pp 最大流問題 pp.85-94
AMPLについて 2011年12月2日(金) 経営システム工学科 森戸 晋.
第7回 条件による繰り返し.
メタ解法設計者のための Python超入門
トーリックイデアルの グレブナ基底を求める アルゴリズム – F4およびF5 –
2. 関係 五島 正裕.
2. 関係 五島 正裕.
スケジューリング最適化システム OptSeq II Pythonモジュールの使い方 補助資料:OptSeq II によるスケジューリング入門 トライアルバージョン
スケジューリング最適化システム OptSeq II 補助資料:OptSeq II によるスケジューリング入門 トライアルバージョン
連続領域におけるファジィ制約充足問題の 反復改善アルゴリズムによる解法 Solving by heuristic repair Algorithm of the Fuzzy Constraint Satisfaction Problems with Continuous Domains 北海道大学.
第4回 線形計画 2000年11月 第4回 線形計画.
部分的最小二乗回帰 Partial Least Squares Regression PLS
ナップサック問題 クマさん人形をめぐる熱いドラマの結末.
情報経済システム論:第13回 担当教員 黒田敏史 2019/5/7 情報経済システム論.
配送計画最適化システム WebMETROのご紹介
制約最適化モジュール SCOP   東京海洋大学 久保幹雄.
サポートベクターマシン Support Vector Machine SVM
データ構造とアルゴリズム (第5回) 静岡大学工学部 安藤和敏
サプライ・チェイン最適化について 研究者・実務家が知っておくべきこと
サプライ・チェイン最適化における モデリングについて
Solver for COnstraint Programming:スコープ Pythonモジュールの使い方
情報工学概論 (アルゴリズムとデータ構造)
平面走査法を使った 一般線分の 交点列挙アルゴリズム
分枝カット法に基づいた線形符号の復号法に関する一考察
情報処理Ⅱ 第7回 2004年11月16日(火).
コストのついたグラフの探索 分枝限定法 A*アルゴリズム.
参考:大きい要素の処理.
情報処理Ⅱ 第2回 2004年10月12日(火).
東京工業大学情報理工学研究科 小島政和 第1回横幹連合コンファレンス 2005年11月25,26日 JA 長野県ビル
Presentation transcript:

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

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

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

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

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 での実装(2)

Python での実装(3) 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 での実装(4) 中略 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 での実装(5) 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) # 描画

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

上界と下界の変化

Optimize a model with Rows, Columns and 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 e+01, best bound e+01, gap 0.0% Opt.value= 強い定式化での結果

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

巡回セールスマン問題(TS P) すべての点をちょうど1回通る最短巡 回路 切除平面法で 85,900 点の実際問題(対称 TSP)の最適解

Miller-Tucker-Zemlin の定式化

AMPL での実装(1) 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];

AMPL での実装(2) 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]; # 持ち上げ上界制約

上界と下界の変化 (80点, Euclid TSP ) 強化した式 でないと... 1日まわして Out of Memory!

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

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

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

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

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

弱い定式化 標準定式化 施設配置定式化 強い定式化 = 線形計画緩和が 整数多面体と一致 変数の数 定式化のサイズと強さの比較 制約の数 変数の数 制約の数 (S,l) 不等式 切除平面 (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 ,彩色数上限 K max =10 ランダムグラフ G(n,p=0.5) 彩色数 8 が最適値

定式化 弱い定式化

Python での実装(1) 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()

Python での実装(2) 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 ,彩色数上限 K max =10 Optimize a model with 3820 Rows, 410 Columns and NonZeros Explored nodes ( simplex iterations) in seconds

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

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

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

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

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

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

制約追加メソッド addConstr の引 数 l hs (左辺):線形制約オブジェクト か定数 sense (制約の向き): GRB.LESS_EQUAL (≦; “ ” でも 可) rhs (右辺):線形制約オブジェクトか 定数 name (名前)

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

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

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

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

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

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