Gurobiの使い方 - LOG OPT HOME

Pythonを用いたGurobi入門
久保 幹雄
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ページ解説
複合型
– リスト:任意の要素から成る順序型
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]:
print i*2
=> 2,4,10,12
for key in D:
print key,D[key]
=>キーと値を出力
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”
変数オブジェクト
Pythonでの実装(2)
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での実装(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:
edge.append((i,j))
return edge
#x[i,j]の解が1なら
#リストedgeに(i,j)を追加
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)
#描画
弱い定式化での結果
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
上界と下界の変化
18
16
Obj. Func. Value
14
12
10
8
6
4
2
0
0
10
20
30
40
CPU
50
60
70
強い定式化での結果
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を用いない強い定式化が望ましい.
• この程度の式なら,必要な式のみを切除
平面として追加するような小細工は必要な
し.
(ただし,式の数が増え退化するので,大
規模なLPを高速に解けるソルバーが前
提)
巡回セールスマン問題(TSP)
• すべての点をちょうど1回通る最短巡回路
• 切除平面法で85,900点の実際問題(対称T
SP)の最適解
Miller-Tucker-Zemlinの定式化
上界と下界の変化
(80点,Euclid TSP)
45
40
35
強化した式
でないと...
1日まわして
Out of Memory!
Obj. Func. Value
30
25
20
15
10
5
0
0
50
100
150
200
CPU
250
300
350
400
結果
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-1)
在庫量(t)
期
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 t
定式化のサイズと強さの比較
標準定式化
変数の数
制約の数
O(n)
施設配置定式化
O(n)
弱い定式化
変数の数 O(n
制約の数
2
追加した
制約の数
O(2n )
O(n )
(S,l)不等式
)
強い定式化
=線形計画緩和が
整数多面体と一致
切除平面
n: 期数
強い定式化
2
上界と下界の変化(標準定式
化)
132000
130000
128000
Obj. Func. Value
126000
124000
122000
120000
118000
116000
114000
112000
0
200
400
600
800
1000
1200
1400
1600
1800
2000
CPU
1800秒で最適解;これ以上大きな問題例は無理!
上界と下界の変化(施設配置定式化)
114050
114000
113950
Obj. Func. Val.
113900
113850
113800
113750
113700
113650
113600
0
5
10
15
20
25
30
35
CPU
40秒で最適解;T=100でも大丈夫!
40
45
知見
• 従来では難問と言われてきたロットサイズ
決定問題でもある程度までは大丈夫
• 緩和固定法(Federgruen, Meissner,Tzur
“Progressive Interval Heuristics for MultiItem Capacitated Lot-Sizing Problems”
2007)というMIPベースのメタ解法を作って
みたが,同時間通常の探索を行う「打ち切
り分枝限定法」の方が良い!
k-center問題
• min-max型の施設配置問題
• 100顧客,10施設(顧客上から選択)
• Euclid距離,x,y座標は一様ランダム
定式化
上界と下界の変化
1.2
Obj. Fun. Value
1
0.8
0.6
0.4
0.2
0
0
50
100
150
200
CPU Time
250
300
350
400
知見
• Min-max型の目的関数はMIPソルバーでは解き
にくい(双対ギャップが大きい;上界も下界も悪い
ので,途中で止めても悪い解!).
• 例: Job shopスケジューリングの最大完了時刻
(メイクスパン)最小化
• 制約計画としてモデル化して,上限を制約として
扱うとうまくいく場合もある.
グラフ彩色問題
•
•
•
•
解の対称性
点数 n=40,彩色数上限 Kmax=10
ランダムグラフ G(n,p=0.5)
彩色数 8 が最適値
定式化
弱い定式化
上界と下界の変化(原定式化)
点数 n=40,彩色数上限 Kmax=10
12
Obj. Func. Value
10
8
6
4
2
0
0
200
400
600
800
1000
1200
CPU Time
Optimize a model with 3820 Rows, 410 Columns and 11740 NonZeros
Explored 17149 nodes (3425130 simplex iterations) in 1321.63 seconds
1400
定式化の改良
対称性の除去(番号の小さい方の色を優先して使う!)
特殊順序集合(SOS: Special Ordered Set) Type 1 (いずれか
1つの変数が正になることの宣言)の追加
model.addSOS(1,変数リスト)
上界と下界の変化(対称性除
去)
12
Obj. Func. Value
10
8
6
4
2
0
0
50
100
150
200
250
300
350
CPU
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秒
400
450
上界と下界の変化(+SOS)
12
Obj. Func. Val.
10
8
6
4
2
0
0
5
10
15
20
CPU
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秒
25
知見
• 解の対称性のある問題は,下界の改善がしにく
いので,分枝限定法では解きにくい(モダンなソ
ルバーは群論などを用いた工夫を入れてはいる
が,あまり機能しない問題もある)
• 対称性を除く工夫を入れると多少は改善
• 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の引数
• lhs(左辺):線形制約オブジェクトか定数
• sense(制約の向き):GRB.LESS_EQUAL
(≦;“<”でも可), GRB.EQUAL(=; “=”で
も可), GRB.GREATER_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)の相対誤差(終了
条件);規定値は 10-4
setParam(“MIPGap”,0.0)
model.Params.MIPGap=0.0
よく使うパラメータ
• TimeLimit:計算時間上限(秒)
• MIPGap:相対誤差上限(規定値は 10-4)
• 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から13まで;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:余裕変数の値