Transcript Power Point - LOG OPT HOME

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ページ解説
複合型
– リスト：任意の要素から成る順序型
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]
=>キーと値を出力
How to solve hard optimization
problems quickly by Python
Metaheuristics
MIP solver Gurobi
Python Source Codes
Developed by:
Me & J. Pedro Pedroso
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={}
キー
“Hanako”,
(1,2)
写像
値
“127cm”
変数オブジェクト
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:
edge.append((i,j))
return edge
#x[i,j]の解が1なら
#リストedgeに(i,j)を追加
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
上界と下界の変化
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を用いない強い定式化が望ましい．
• この程度の式なら，必要な式のみを切除
平面として追加するような小細工は必要な
し．
（ただし，式の数が増え退化するので，大
規模なＬＰを高速に解けるソルバーが前
提）
k-center問題
• min-max型の施設配置問題
• 施設は顧客上から選択
• Euclid距離，x,y座標は一様ランダム
Which is the solution of the k-center problem?
定式化
最大距離を最小化
上界と下界の変化
1.2
Obj. Fun. Value
1
0.8
0.6
0.4
0.2
0
0
50
100
150
200
250
300
350
CPU Time
100顧客，10施設 （k-medianは200施設，20顧客で3秒）
400
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-被覆問題を用いた二分探索だとある程度minmax型の目的関数を回避できる．
• min-max型の問題は，制約計画としてモデル化し
て，上限を制約として扱うと良い．
巡回セールスマン問題（ＴＳＰ）
• すべての点をちょうど１回通る最短巡回路
• 切除平面法で85,900点の実際問題（対称Ｔ
ＳＰ）の最適解
Miller-Tucker-Zemlinの定式化
上界と下界の変化
（８０点，Euclid TSP）
45
40
35
強化した式
でないと．．．
１日まわして
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
単一フロー定式化
多少強化
対称な問題（DFJ定式化）
部分巡回路除去制約を追加
線形緩和問題
部分巡回路除去制約を追加
整数解
連結成分を求めて切除平面を追加
部分巡回路除去制約を追加
最適解
知見
• 式を持ち上げ操作などで強化すると，高速化さ
れ，大きな問題例が解けるようになる．
（そのためには多面体論の知識が多少必要なの
で，専門家に相談．ちょっと式をいじるだけで劇
的に改善する場合もある！）
• 他にも単品種フロー，多品種フロー定式化など
がある．（下界が強い方が良いとは限らない！
LPのサイズと解きやすさ（退化の具合）を考慮し
て定式化を選択！）
• 対称な問題なら分枝カット法が良い．
多品目ロットサイズ決定問題
• 段取り費用と在庫費用のトレードオフを最
適化する多期間生産計画
• 多品目で共通の資源を使う容量制約付き
問題は，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ベースのメタ解法を作って
みたが，同時間通常の探索を行う「打ち切
り分枝限定法」の方が良い！
ビンパッキング問題
6
6
7
5
5
5
4
4
4
8
7
5
4
2
2
2
2
3
5
8
4
3
4
5
｝
9
アイテム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)
master.update()
Pythonによる記述（列生成）
while 1:
knapsack.update()
iter += 1
knapsack.optimize()
relax = master.relax()
if knapsack.ObjVal < 1+EPS:
relax.optimize()
break
pi = [c.Pi for c in relax.getConstrs()]
col = Column()
knapsack = Model("KP")
for i in range(m):
knapsack.ModelSense=-1
if t[K][i] > 0:
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,
ame="width")
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()
列生成法の応用事例
（タンカースケジューリング問題）
• タンカーのパーティションを考慮した複雑な
実際問題
• 海技研からの委託
• 可能なパターンの列挙，集合被覆問題に
よる選択を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))
for k in range(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
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 （いずれか
１つの変数が正になることの宣言）の追加
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
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
スカラー化が見逃す非劣解
第一目的関数
二目的巡回セールスマン問題に対する制約スライド法の結果
制約スライド法と理想点法の比較
Gurobi Objects
GRBError
addVar
Variable
Model
addSOS
addConstr
Column
Constraint
SOS
Callbacks
LinExpr
QuadExpr
Model Object
= Model(“name of model”）
Methods
• AddVar
• AddConstr
• optimize
• update
• getVars
• relax
...
Attributes
•
•
•
•
•
•
...
Params
ObjVal
ModelSense
Runtime
Status
SolCount
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
SCOP Objects
addVariable(s)
Model
addConstraint
Variable
Linear
Quadratic
Alldiff
詳細については http://www.logopt.com/scop/
Model Object
= Model(“name of model”）
Methods
• AddVariable(s)
• AddConstraint
• optimize
Attributes
Params
• variables
• constraints
•
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
OptSeq Objects
addActivity
Attribute
addMode
Model
addResource
addTemporal
Mode
Resource
Temporal
詳細については http://www.logopt.com/OptSeq/
Model Object
= Model(“name of model”）
Methods
•
•
•
•
•
AddActivity
AddResource
AddTemporal
optimize
write
Attributes
•
•
•
•
•
Params
activities
modes
resources
temporals
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