Transcript 講演資料
最適化ワークショップ 整数計画法の発展と応用
整数計画法チュートリアル:モデリングと解法
藤江哲也
兵庫県立大学大学院経営研究科
2013年5月22日(水) 九州大学伊都キャンパス
内容
整数計画問題とは
整数計画問題の例(ナップサック問題、分割問題)
整数計画法のShort History,分枝限定法
整数計画による定式化について
おわりに
2
線形計画問題
A社では,テーブルとチェアを製造販売している.それぞれ製造工程が2つ
あり,1個当たり所要時間、利益、および製造工程の使用可能時間が次のよ
うに与えられている.
テーブル
チェア
使用可能時間
工程1
2時間
2時間
7時間
工程2
3時間
5時間
14時間
利益
4千円
5千円
テーブルとチェアはすべて売れるものとした場合,利益を最大にするにはテ
ーブルとチェアをそれぞれ1日当たり何個製造すればよいか.
3
線形計画問題
x1
x2
テーブル
チェア
使用可能時間
工程1
2時間
2時間
7時間
工程2
3時間
5時間
14時間
利益
4千円
5千円
最大化 z 4 x1 5 x2
条件
2 x1 2 x2 7
3 x1 5 x2 14
x1 , x2 0
x2
最大化
2 x1 2 x2 7
3 x1 5 x2 14
z 4 x1 5 x2 15.75
z 4 x1 5 x2 10
最適解 (x1, x2) = (7/4, 7/4)
2
z = 63/4 = 15.75
z 4 x1 5 x2 5
z 4 x1 5 x2 0
1
O
1
2
3
x1
4
整数計画問題(整数線形計画問題)
x1
x2
テーブル
チェア
使用可能時間
工程1
2時間
2時間
7時間
工程2
3時間
5時間
14時間
利益
4千円
5千円
最大化 z 4 x1 5 x2
条件
x2
最大化
2 x1 2 x2 7
3 x1 5 x2 14
x1 , x2 0
x1 , x2 : 整数
最適解 (x1, x2) = (1, 2)
2
z = 14
1
O
1
2
3
x1
5
混合整数計画問題
MIP (Mixed Integer Programming)
整数条件:一部またはすべての変数
最大化 z 4 x1 5 x2
条件
x2
最大化
2 x1 2 x2 7
3 x1 5 x2 14
x1 , x2 0
x1 : 整数
最適解 (x1, x2) = (2, 3/2)
2
z = 31/2=15.5
1
O
1
2
3
x1
6
混合整数計画問題(MIP Problems)
線形計画問題+「変数の整数条件」
n
最小化
c x
j 1
n
条件
j
または
j
a x
j 1
ij
j
bi
xj 0
x j : 整数
(i 1,, m)
( j 1, , n)
( j 1, , p )
最小化 cT x
条件
Ax b
x0
x Z p R n p
( p n)
7
線形計画と整数計画
P(多項式時間)
線形計画はクラス
代表的解法
単体法(シンプレックス法)
内点法
個人用PCでも、数万変数、数万制約以上の問題を解くこ
とができる
NP困難
整数計画は
8
混合整数計画問題(MIP Problems)
0-1変数(バイナリ変数)
『 𝑥𝑗 = 0 または 1 』
「 𝑥𝑗 = 0 または 1 」 ⟺ 「0 ≤ 𝑥𝑗 ≤ 1, 𝑥𝑗 ∈ 𝑍」より、
0-1変数も 「整数変数+線形不等式」 で記述できる
しかし、 0-1変数は整数変数と区別されるのが一般的
高い表現能力
0-1の特性に基づくアルゴリズム開発
9
例)ナップサック問題
アイテム j
1
2
3
4
5
6
価値 pj
10
3
16
2
10
3
重さ wj
3
4
7
2
4
5
重さの合計が容量以下という制約の下で,価
値の合計が最大になる組み合わせを求めよ
重さの合計 = 4 + 7 + 2 = 13 <= 14
価値の合計 = 3 + 16 + 2 = 21
ナップサックの容量
c = 14
10
例)ナップサック問題
アイテム j
1
2
3
4
5
6
価値 pj
10
3
16
2
10
3
重さ wj
3
4
7
2
4
5
1 アイテム j を選ぶとき
xj
0 アイテム j を選ばないとき
とすると
最大化 10 x1 3 x2 16 x3 2 x4 10 x5 3 x6
条件
3 x1 4 x2 7 x3 2 x4 4 x5 5 x6 14
x1 , x2 , , x6 0 または 1
x1=1 のとき 3
x1=0 のとき 0
11
例)ナップサック問題
一般に
最大化
p x
w x
jJ
条件
jJ
j
j
j
j
c
x j 0 または 1 ( j J )
x j 0, 1 ( j J )
と書くことも多い
12
例)分割問題
(1,2,5,7,8,9) を,合計が等しくなるように2つの集合に分割
できるか?
一般に
(𝑎1 , 𝑎2 , ⋯ , 𝑎𝑛 ) を,合計が等しくなるように2つの集合に分
割できるか?
𝐴=
𝑛
𝑖=1 𝑎𝑖
とするとき,
Find 𝐼 ⊆ {1,2, ⋯ , 𝑛} s.t.
𝑖∈𝐼 𝑎𝑖
=𝐴 2
13
例)分割問題
1
xj
0
j を選ぶとき ( j I )
とする
j を選ばないとき ( j I )
方法1
最大化 0 x1
n
条件
a x
j 1
j
j
A/ 2
x j {0,1} ( j 1, , n)
方法2
最小化 s
n
最大化
a x
j 1
n
条件
j
j
ajxj A/ 2
j 1
x j {0,1} ( j 1, , n)
n
条件
a x
j 1
j
j
s A/ 2
x j {0,1} ( j 1, , n)
s0
14
例)分割問題
(1,2,5,7,8,9) を,合計が等しくなるように3つの集合に分割
することができるか?
⇔ スケジューリング問題 (P||Cmax)
機械3
9
機械2
機械1
5
1
2
8
7
13
15
例)分割問題
最小化 s
n
a x
条件
j 1
m
j ij
s
n
最小化 max a j xij
j
j 1
(i 1, , m)
m
xij 1
条件
( j 1, , n)
i 1
x j {0,1} (i 1, , m; j 1, , n)
1
(j )
8
機械2
機械1
ij
x j {0,1} (i, j )
i 1
機械3
x
5
1
2
7
9
s=12
16
例)分割問題
最大化 t
n
a x
条件
j 1
m
j ij
t
n
最大化 min a j xij
j
j 1
(i 1, , m)
m
xij 1
条件
( j 1, , n)
i 1
x j {0,1} (i 1, , m; j 1, , n)
機械2
機械1
1
ij
1
(j )
x j {0,1} (i, j )
i 1
機械3
x
9
2
8
5
7
t=10
12
17
例)分割問題
最小化 s t
n
条件
a x
j ij
j 1
n
a x
j 1
m
j ij
x
i 1
ij
s
n
n
最小化 max a j xij min a j xij
j
j
j 1
j 1
(i 1, , m)
m
t
1
条件
(i 1, , m)
x
i 1
ij
1
(j )
x j {0,1} (i, j )
( j 1, , n)
x j {0,1} (i 1, , m; j 1, , n)
機械3
機械2
機械1
1
9
8
2
5
7
t=10
s=12
18
例)分割問題
実行可能解を求める問題(実行可能性判定問題)
⇒ 最適化問題として考えることができる(cf. SATとMaxSat)
※最適化問題 → 実行可能性判定
𝑧∗ =
𝑐𝑗 𝑥𝑗∗ の解が既知 → 𝑧 ∗ <
𝑐𝑗 𝑥𝑗 を満たす解が存在するか?
実問題では,制約をすべてを満たすことが困難/無理な場合
が多く見られる
⇒ ハード制約とソフト制約に分類
⇒ ソフト制約をなるべく満たす ⇒ 最適化問題
分割問題その他の目的関数(評価基準)?
機械ごとに重みをつける?
分散最小化?
19
2次関数の線形化
0-1変数の場合、2次関数を線形化することができる
・ 𝑥𝑗2 = 𝑥𝑗
𝑦𝑖𝑗 = 𝑥𝑖 𝑥𝑗 ,
・ 𝑥 , 𝑥 ∈ {0,1} ↔
𝑖 𝑗
1 − 𝑥𝑖 − 𝑥𝑗 + 𝑦𝑖𝑗 ≥ 0,
𝑥𝑖 − 𝑦𝑖𝑗 ≥ 0, 𝑥𝑗 − 𝑦𝑖𝑗 ≥ 0,
𝑦𝑖𝑗 ≥ 0, 𝑥𝑖 , 𝑥𝑗 ∈ {0,1}
1
yij
1
1
xi
同様に多項式も線形化できる
0-1変数を使うと,一般の非線形関数 𝑦 = 𝑓(𝑥) の近似がで
きる(次ページ)
20
xj
非線形関数 𝑦 = 𝑓 𝑥 の区分線形近似
𝑥=
𝜆𝑘 𝑥𝑘 , 𝑦 =
𝜆𝑘 𝑓 𝑥𝑘 ,
𝜆𝑘 = 1, 𝜆𝑘 ≥ 0 ∀𝑘 ,
𝜆1 ≤ 𝑧1 ,
𝜆𝑘 ≤ 𝑧𝑘−1 + 𝑧𝑘 2 ≤ 𝑘 ≤ 𝐾 − 1 ,
𝑧𝑘 ≤ 𝑧𝐾−1 ,
𝑧𝑘 = 1, 𝑧𝑘 ∈ {0,1} ∀𝑘
(𝑥𝐾 , 𝑓 𝑥𝐾 )
(𝑥2 , 𝑓 𝑥2 )
y
(𝑥1 , 𝑓 𝑥1 )
x
21
非凸集合の表現
離接(disjunctive)制約
𝑥1 + 2𝑥2 ≤ 2 または 2𝑥1 + 𝑥2 ≤ 2
x1 2 x2 2 M (1 y )
2 x1 x2 2 My
y {0,1}
big-M
(十分大きい数)
x2
2
1
O
1
2
x1
22
その他0-1変数を用いた表現
0-1変数 : 「採用/不採用」 「Yes/No」 「true/false」 ...
例1) 𝑎 ∨ 𝑏 ∨ c
𝑎+ 1−𝑏 +𝑐 ≥1
例2) a, b のうち高々1つを採用
𝑎+𝑏 ≤1
例3) 𝑥1 , … , 𝑥𝑛 のうち高々𝑘個を採用
𝑛
𝑗=1 𝑥𝑗
≤𝑘
例4) aを採用 → bを採用
𝑎≤𝑏
23
その他0-1変数を用いた表現
𝑧
固定費付き
𝑧 = 𝑑𝑥 + 𝑓𝑦, 0 ≤ 𝑥 ≤ 𝑀𝑦, 𝑦 = 0 または 1
𝑧 = 𝑑𝑥 + 𝑓
big-M
𝑥
施設配置
𝑥𝑗 ≤ 𝑀𝑦, 𝑥𝑗 ≥ 0 ∀𝑗 , 𝑦 = 0 または 1
big-M
open
or
not
𝑥𝑗
24
整数計画問題による定式化
“文章題”による例題では適用範囲が狭い印象を与える
しかし、実際には幅広い適用可能性がある
0-1変数
非線形、非凸の問題をカバー
実行可能解を求める問題(実行可能性判定)を最適化問題として捉
える
定式化はいろいろ
問題の捉え方による表現の違い
同じ問題でも...(→ 後述)
25
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、...
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
組合せ最適化問題のように、問題固有の構造を持たない
組合せ最適化問題:実行可能解を求めるのは容易な場合が多い
整数計画問題
:実行可能性判定も一般にNP完全
問題によって
は現在も必要
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、1990年代以降、整数計画ソルバーが目覚ましく進展
④ それに伴い、整数計画が利用しやすくなっている
26
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、...
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
組合せ最適化問題のように、問題固有の構造を持たない
組合せ最適化問題:実行可能解を求めるのは容易な場合が多い
整数計画問題
:実行可能解を求めることも一般にNP完全
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、1990年代以降、整数計画ソルバーが目覚ましく進展
④ それに伴い、整数計画が利用しやすくなっている
27
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、...
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
組合せ最適化問題のように、問題固有の構造を持たない
組合せ最適化問題:実行可能解を求めるのは容易な場合が多い
整数計画問題
:実行可能解を求めることも一般にNP完全
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、1990年代以降、整数計画ソルバーが目覚ましく進展
④ それに伴い、整数計画が利用しやすくなっている
28
② Short History
1958年 R. E. Gomory
"Outline of an Algorithm for Integer Solutions to Linear Programs",
Bulletin of the American Mathematical Society
全整数計画問題
小数カットに基づく切除平面法
(cutting plane algorithm)
有限回の繰り返しで最適解が得られる
最小化
LP解
LP解
LP解
=最適解!
29
② Short History
1960年 R. E. Gomory
“An algorithm for the mixed integer problem”, Technical Report RM-2597,
The Rand Corporation
混合整数計画問題 (MIP)
カットの提案 (Gomory Mixed Integer cut : GMI cut)
1960年 A. Land and A. Doig
"An Automatic Method for Solving Discrete Programming Problems", Econometorica
など
分枝限定法の提案
30
分枝限定法
緩和法の原理
元問題
(MIP ) 最小化 c T x
条件
Ax b
x0
x Z p R n p
LP緩和問題
(LP ) 最小化 c T x
条件
Ax b
x0
x Rn
31
分枝限定法
(MIP)の実行可能領域
(LP)の実行可能領域
[定理(緩和法の原理)]
(R1) (MIP)の最適値
(LP)の最適値
(R2) (LP)の最適解 x* が整数条件を満たす ⇒ x* は(MIP)の最適解
(R3) (LP)が実行不能 ⇒ (MIP)も実行不能
目的関数値
(最小化)
(MIP) の実行可能解
暫定値(上界値)
下げたい
(MIP) の最適値
上げたい
(LP) の最適値
下界値
32
分枝限定法
(LP)を解く.(LP)の最適解 x* が整数条件を満たせば終了(R2より).
さもなければ変数 xt を選び, xt xt* に固定した問題 (MIP1) と
xt xt* に固定した問題 (MIP2) を生成する(分枝操作)
(MIP ) 最小化 c T x
条件
Ax b
x0
x Z p R n p
例) x3 に関する分枝操作
MIP
x3 2
x3 3
(MIP1 ) 最小化 cT x
条件
Ax b
x0
x3 2
x Z p R n p
(MIP2 ) 最小化 cT x
MIP1
MIP2
条件
Ax b
x0
x3 3
x Z p R n p
33
分枝限定法
限定操作
子問題 (MIPk) に対する限定操作
(MIPk)の緩和問題(LPk)を解く
1
(LPk)が実行不能 ⇒ (MIPk) は実行不能 (R3より)
2
(LPk) の最適解 xk ,最適値 zk とする
2-1
xk が整数 ⇒
xk は (MIPk) の最適解 (R2より)
よって,zk < z* (暫定値) ならば,z* = zk と更新する
2-2
3
zk ≥ z* (暫定値) ならば, (MIPk) の中に最適解はない(R1より)
1 でも 2 でもなければ分枝操作を行う
34
分枝限定法
例題
(MIP) 最小化 2 x1 3 x2
条件
x2
最小化
2 x1 x2 10
3 x1 6 x2 40
x1 0, x2 0
x1 , x2 : 整数
(LP) 最小化 2 x1 3x2
条件
2 x1 x2 10
3x1 6 x2 40
x1 0, x2 0
O
x2
x1
20 50
,
,
9 9
190
𝑧=−
= −21.1
9
𝑥, 𝑦 =
下界値
O
x1
35
分枝限定法
分枝限定木(列挙木)
z1 = -21
(2,17/3)
P
x1 ≤ 2
P2
x2 ≥ 6
P3
P4
z*=-19
z6 = -20
(1,6)
z*=-20
x2 ≤ 6
P7
z2 = -18
(3,4)
z4 = -62/3
(4/3,6)
x1 ≤ 1
z5 = -41/2
(1,37/6)
z*=+∞
x1 ≥ 3
P1
x2 ≤ 5
z3 = -19
(2,5)
z = -190/9
(20/9,50/9)
最適解 (1,6)
最適値 -20
x1 ≥ 2
P5
P6
実行不能
x2 ≥ 7
P8
実行不能
36
分枝限定法
目的関数値
暫定値(上界値)
下界値
時間
37
分枝限定法
分枝前の
MIP
x3 2
MIP1
z1 z
x3 3
MIP2
z2 z
下界値 z
分枝後の
下界値 min z1, z2 z
一般に・・・
(グローバルな )下界値
minノード i の下界値
i
38
分枝限定法
近年は分枝カット法(分枝限定法+カット)が主流
分枝限定法/分枝カット法の実装にあたって
子問題の選択 (Node selection)
下界値優先 (Best Bound) :下界値の良い(小さい)子問題を優先
評価値優先 (Best Estimate):評価値の良い子問題を優先
深さ優先 (Depth First)
:分枝されたばかりの子問題を優先する
分枝変数の選択 (Variable selection)
前処理 (Preprocessing)による定式化の強化
カット (Cutting plane)
ヒューリスティック (Heuristic)
・・・
39
② Short History
ともに線形計画をベースとする解法
実用的には分枝限定法
切除平面法 : 最適解への収束が遅い、数値誤差がたまる、解法途中で実行可能
解が得られない
組合せ最適化問題ごとに解法を開発
※巡回セールスマン問題
・1954年 G. Dantzig, R. Fulkerson and S. M. Johnson,
“Solution of a Large-scale Traveling-salesman", Operations Research
・1990年 M. Padberg and G. Rinaldi,
“A Branch-and-cut algorithm for the Resolution of Large-scale Symmetric Traveling
Salesman Problems", SIAM Review
・1997年 D. L. Applegate, R. E. Bixby, V. Chvátal and W. J. Cook,
“The Traveling Salesman Problem: A Computational Study”, Princeton University Press
40
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、...
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
組合せ最適化問題のように、問題固有の構造を持たない
組合せ最適化問題:実行可能解を求めるのは容易な場合が多い
整数計画問題
:実行可能解を求めることも一般にNP完全
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、1990年代以降、整数計画ソルバーが目覚ましく進展
④ それに伴い、整数計画が利用しやすくなっている
41
③ 1990年代に整数計画ソルバーが急速に進展
計算機パワーの増大
高速でロバストな線形計画ソルバーの開発
分枝カット法=分枝限定法+切除平面法
1990年 M. Padberg and G. Rinaldi (巡回セールスマン問題)
Gomoryカット(GMI cut)の適用
1996年 E. Balas, S. Ceria, G. Cornuejols and N. Natraj,
“Gomory Cuts Revisited", Operations Research Letters
42
③ ソルバーの進展(分枝限定法の効率化)
子問題の選択
分枝変数の選択
前処理による定式化の強化
カット
ヒューリスティクス
Constraint Programming との融合
対称性
優越テスト
並列化(マルチスレッド)
※ロバストな線形計画ソルバーが必要
43
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、...
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
組合せ最適化問題のように、問題固有の構造を持たない
組合せ最適化問題:実行可能解を求めるのは容易な場合が多い
整数計画問題
:実行可能解を求めることも一般にNP完全
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、1990年代以降、整数計画ソルバーが目覚ましく進展
④ それに伴い、整数計画が利用しやすくなっている
44
④ 整数計画ソルバーを使用する
商用
FICO Xpress (Fair Isaac Corporation)
Gurobi Optimizer (Gurobi Optimization)
IBM ILOG CPLEX (IBM)
LINDO (LINDO Systems)
NUOPT (数理システム)
SOPT (Saitech, Inc.)
Excel Solver (Frontline Systems)
・・・・・
フリー
COIN/CBC
GNU GLPK
lp_solve
SCIP
・・・・・
: H. Mittelmannによる
ベンチマーク比較
45
定式化について
定式化とは?
⇒ 実行可能解集合を正しく表現する「不等式系+整数条件」
46
定式化について
定式化は複数ありうる
変数の定義の仕方によっても(大きく)異なる
(復習)現在の整数計画ソルバー : 線形計画ベースの分枝
カット法
⇒ 緩和の強い定式化が望ましい
現在では、緩和を強化する機能を備えた整数計画ソルバー
も多い
前処理(Presolve, Preprocessing)
緩和の強化に加え、不要な変数・制約式の除去なども行う
47
定式化について
とは言え、緩和の強い定式化が望ましい
※その意味で、big-Mを含む定式化は慎重に
次のスライド
また、変数・制約式が少ない定式化が望ましい
どちらを優先的に考えるかはケースバイケース
48
Big-Mを含む問題例
x2
最大化 𝑧 = 2𝑥1 + 3𝑥2
条 件 𝑥1 + 2𝑥2 ≤ 2 または 2𝑥1 + 𝑥2 ≤ 2
𝑥1 , 𝑥2 ≥ 0
2
最大化 𝑧 = 2𝑥1 + 3𝑥2
1
条 件 𝑥1 + 2𝑥2 − 2 ≤ 𝑀(1 − 𝑦)
2𝑥1 + 𝑥2 − 2 ≤ 𝑀𝑦
𝑥1 , 𝑥2 ≥ 0, 𝑦 ∈ 0,1
O
1
2
x1
MIPの最適解 𝑥1 , 𝑥2 , 𝑦 = 0, 2, 1 , 𝑧 = −6
𝑀 = 10000の場合
LPの最適解 𝑥1 , 𝑥2 , 𝑦 = 0, 3334.67, 0.33 , 𝑧 = −10004
49
定式化の例1:数独
7
5
9
6
5
3
3
4
3
6
9
1
3
1
3
2
1
9
5
8
6
4
3
4
5
9
7
6
3
1
2
8
6
3
8
4
1
2
5
7
9
1
7
3
6
8
4
9
5
2
5
9
4
1
2
7
3
8
6
7
8
6
2
3
9
5
4
1
7
5
2
4
6
8
3
1
7
9
5
3
8
7
5
4
9
2
6
1
9
1
5
2
7
6
8
3
4
2
5
1
9
9
2
7
2
4
7
1
8
9
8
4
7
7
6
4
• あいているマスに 1-9 までのどれかの数字を入れる。
• 縦・横の各列及び、太線で囲まれた3×3のブロックに同じ数字が入ってはいけない。
Wikipediaの「数独」より
50
数独の定式化
方法1
xij : (i, j)マスに入る数字(つま
り1 xij 9)
方法2
1 (i, j )マスに入る数字がk のとき
xijk
0 そうでないとき
T. Koch, "Rapid Mathematical Programming or How to Solve Sudoku Puzzles in a Few Seconds",
Operations Research Proceedings 2005
51
定式化:方法1
xij : (i, j)マスに入る数字(つま
り1 xij 9)
最小化 arbitrarily
条件
xij xi 1
xij xkj 1
x3( r 1) i ,3( c 1) j x3( r 1) k ,3( c 1) 1
1 xij 9
xij k
xij : 整数
i, j, 1,, n; j
i, j, k 1,, n; i k
i, j, k , 1,2,3; 3i j 3k
(i, j 1,, n)
初期配置で(i, j )マスが k のとき
i, j 1, n
n=9
52
定式化:方法1
注1
xij xi x x
b x 8b
xij xi 1 b x 8b
b b 1
b , b 0 または1
注2
constraint programmingの分野では
all_different( xi1,, xin )
という表現をする
53
定式化:方法2
1 (i, j )マスに入る数字がk のとき
xijk
0 そうでないとき
最小化 arbitrarily
n
条件
x
1
i, j 1,, n
x
1
i, k 1,, n
x
1
j, k 1,, n
ijk
k 1
n
ijk
j 1
n
ijk
i 1
3r
3c
x
ijk
i 3( r 1) 1 j 3( c 1) 1
xijk 1
xijk 0 または1
1
r , c 1,2,3; k 1, n
初期配置で(i, j )マスが k のとき
i, j, k 1, n
54
定式化の比較
定式化1
変数 3970(一般整数変数 2025、0-1変数 1944)
制約 4860+固定制約数
計算時間 41.5秒(Intel Core2 Duo [email protected], メモリ2.0GB,
CPLEX 12.5)
※Kochの論文によると、CPLEX9.03では6時間経っても解けなかった
定式化2
変数 729(0-1変数 729)
制約 346+固定制約数
計算時間 0.02秒
55
巡回セールスマン問題(TSP)
Traveling Salesman Problem
各都市をちょうど1回ずつ訪問し、
出発地点に戻ってくる閉路(巡回
路)の中で、ウェイト(距離、時間
等)の合計が最小となるものを求
めよ
兵庫県TSP
56
(非対称)巡回セールスマン問題
Asymmetric Traveling Salesman Problem(ATSP)
1
2
都市数:n
6
都市 i から j への移動距離:cij
cij
3
5
4
j
i
cij = cji を仮定しない
cji
1
巡回路
2
6
3
5
1 : (i, j )が巡回路に含まれると
き
xij
0 : 含まれないとき
として定式化
4
57
定式化(Dantzig-Fulkerson-Johnson, 1954)
(ATSP-DFJ)
n
最小化
n
c x
i 1 j 1
n
条件
x
j 1
ij
n
x
i 1
ij
x
i , jS
ij
ij ij
1
(i 1, , n)
1
( j 1, , n)
S 1 ( S V , 1 S , S 1)
xij 0 または1 (i, j 1, , n)
出る辺は1本
i
j
入る辺は1本
部分巡回路
除去制約
58
部分巡回路除去制約
Dantzig-Fulkerson-Johnson (1954) 再掲
xij | S | 1 (S V , 1 S , | S | 1)
i , jS
変数 𝑂 𝑛2
制約 𝑂(2𝑛 )
部分巡回路(subtour)
S {1, 2}
1
2
6
x12 x21 1
S {3, 4, 5, 6}
3
5
x34 x43 x35 x56 x65 3
4
59
部分巡回路除去制約の表現
Miller-Tucker-Zemlin (1960)
(ATSP-MTZ)
u1 1
2 ui n
(i 2,, n)
変数 𝑂 𝑛2
ui u j 1 (n 1)(1 xij )
(i, j 2,, n)
制約 𝑂(𝑛2 )
ui
・・・ノード i を訪れる順番示す変数
変数の数は n 増加したが,
制約式の数が大幅に減少
60
部分巡回路除去制約の表現
Sherali-Driscoll (2002)
yij (n 1) xi1 ui
(i 2,, n)
(ATSP-SD)
y
( j 2,, n)
変数 𝑂 𝑛2
(i, j 2,, n)
制約 𝑂(𝑛2 )
j 1
i 1
ij
1 u j
xij yij (n 2) xij
u j (n 2) xij (n 1)(1 xij )
yij y ji ui (1 x ji )
(i, j 2,, n)
1 (1 x1 j ) (n 3) x j1
u j (n 1) (n 3) x1 j (1 x j1 )
yij ui xij
( j 2,, n)
実際は変数の数も制約式
の数も増えている
・・・ (ATSP-MTZ)の強化
61
部分巡回路除去制約の表現
Gouveia-Pires (1999, 2001)
xij yij
(i, j 2, , n)
(ATSP-GP)
yij y ji 1
(i, j 2, , n)
変数 𝑂 𝑛2
xij y ji 1
(i, j 2, , n)
制約 𝑂(𝑛3 )
xij yki ykj 1
(i, j , k 2, , n) (☆)
yij 0 または 1
(i, j 2, , n)
1 : i が j に先行するとき
yij
0 : 先行しないとき
(ATSP-GP1)
(ATSP-GP2)
(☆)の
強化
(ATSP-GP3)
62
部分巡回路除去制約の表現
Wong (1980)
yijk xij
n
y
(i, j 1, , n; k 2, , n)
(ATSP-FL)
1
(k 2, , n)
変数 𝑂 𝑛2
yik1 0
(k 2, , n)
制約 𝑂(𝑛3 )
i 2
k
1i
n
i 2
n
y
i 1
k
ik
1
(k 2, , n)
k
kj
0
(k 2, , n)
n
y
j 1
n
n
y y
i 1
k
ij
yij 0
k 2
k
ji
0
( j , k 2, , n)
(i, j 1, , n)
yijk : flow量 (k - th unit)
63
定式化の比較(att48:48都市)
Intel Core2 Duo [email protected], メモリ2.0GB, CPLEX 12.2
変数
制約式
初期LP
IP最終解
目的関数値
時間(秒)
上界値
下界値
誤差
時間(秒)
MTZ
2304
2258
8490.11
0.01
10628
10377.77
2.35% >1800
SD
4465
8932
10070.29
1.43
10628
10523.99
0.98% >1800
GP
4418
101710
9954.50
0.27
11042
10136.65
8.20% >1800
GP1
4418
101757
10041.50
0.49
10685
10313.33
3.48% >1800
GP2
4418
101757
10053.50
0.47
10719
10308.47
3.83% >1800
GP3
4418
199047
10053.50
1.75
11072
10204.43
7.84% >1800
112848
108384
10604.00
236.95
10628
10628.00
0.00%
FL
628.00
64
余談
Excel Solverの「all different」を使った解法(att48の例)
65
余談
巡回セールスマン問題は「Concorde」が有名
66
おわりに
まとめ1
整数計画問題:「解けない」から「解ける」へ
MIP:汎用的 → インパクトが大きい
現在のMIPソルバー = LPベースの分枝カット法
⇒ 苦手な問題の種類もある(big-Mを含む問題など)
まとめ2
問題の捉え方によって定式化が異なる
変数の設定方法によって定式化が異なる
⇒「悪い定式化」になる可能性も
しかし、まずは定式化をしてみる
安心のため...
強引な定式化でも、今のソルバーなら解けてしまうかも
67