Transcript 講演資料

最適化ワークショップ 整数計画法の発展と応用
整数計画法チュートリアル：モデリングと解法
藤江哲也
兵庫県立大学大学院経営研究科
2013年5月22日（水） 九州大学伊都キャンパス
内容
 整数計画問題とは
 整数計画問題の例（ナップサック問題、分割問題）
 整数計画法のShort History，分枝限定法
 整数計画による定式化について
 おわりに
2
線形計画問題
A社では，テーブルとチェアを製造販売している．それぞれ製造工程が２つ
あり，１個当たり所要時間、利益、および製造工程の使用可能時間が次のよ
うに与えられている．
テーブル
チェア
使用可能時間
工程１
2時間
2時間
7時間
工程２
3時間
5時間
14時間
利益
4千円
5千円
テーブルとチェアはすべて売れるものとした場合，利益を最大にするにはテ
ーブルとチェアをそれぞれ１日当たり何個製造すればよいか．
3
線形計画問題
x1
x2
テーブル
チェア
使用可能時間
工程１
2時間
2時間
7時間
工程２
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
テーブル
チェア
使用可能時間
工程１
2時間
2時間
7時間
工程２
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
x0
x  Z p  R n p
( p  n)
7
線形計画と整数計画
Ｐ（多項式時間）
 線形計画はクラス
 代表的解法
単体法（シンプレックス法）
内点法
 個人用ＰＣでも、数万変数、数万制約以上の問題を解くこ
とができる
ＮＰ困難
 整数計画は
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
jJ
条件
jJ
j
j
j
j
c
x j  0 または 1 ( j  J )
x j  0, 1 ( j  J )
と書くことも多い
12
例）分割問題
(1,2,5,7,8,9) を，合計が等しくなるように２つの集合に分割
できるか？
一般に
(𝑎1 , 𝑎2 , ⋯ , 𝑎𝑛 ) を，合計が等しくなるように２つの集合に分
割できるか？
𝐴=
𝑛
𝑖=1 𝑎𝑖
とするとき，
Find 𝐼 ⊆ {1,2, ⋯ , 𝑛} s.t.
𝑖∈𝐼 𝑎𝑖
=𝐴 2
13
例）分割問題
1
xj  
0
j を選ぶとき ( j  I )
とする
j を選ばないとき ( j  I )
方法１
最大化 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)
s0
14
例）分割問題
(1,2,5,7,8,9) を，合計が等しくなるように３つの集合に分割
することができるか？
⇔ スケジューリング問題 （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
２次関数の線形化
 0-1変数の場合、２次関数を線形化することができる
・ 𝑥𝑗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」 ．．．
例１） 𝑎 ∨ 𝑏 ∨ c
𝑎+ 1−𝑏 +𝑐 ≥1
例2） a, b のうち高々１つを採用
𝑎+𝑏 ≤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
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、．．．
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
 組合せ最適化問題のように、問題固有の構造を持たない
 組合せ最適化問題：実行可能解を求めるのは容易な場合が多い
整数計画問題
：実行可能性判定も一般にＮＰ完全
問題によって
は現在も必要
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、1990年代以降、整数計画ソルバーが目覚ましく進展
④ それに伴い、整数計画が利用しやすくなっている
26
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、．．．
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
 組合せ最適化問題のように、問題固有の構造を持たない
 組合せ最適化問題：実行可能解を求めるのは容易な場合が多い
整数計画問題
：実行可能解を求めることも一般にＮＰ完全
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、1990年代以降、整数計画ソルバーが目覚ましく進展
④ それに伴い、整数計画が利用しやすくなっている
27
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、．．．
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
 組合せ最適化問題のように、問題固有の構造を持たない
 組合せ最適化問題：実行可能解を求めるのは容易な場合が多い
整数計画問題
：実行可能解を求めることも一般にＮＰ完全
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、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
x0
x  Z p  R n p
LP緩和問題
(LP ) 最小化 c T x
条件
Ax  b
x0
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
x0
x  Z p  R n p
例） x3 に関する分枝操作
MIP
x3  2
x3  3
(MIP1 ) 最小化 cT x
条件
Ax  b
x0
x3  2
x  Z p  R n p
(MIP2 ) 最小化 cT x
MIP1
MIP2
条件
Ax  b
x0
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
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、．．．
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
 組合せ最適化問題のように、問題固有の構造を持たない
 組合せ最適化問題：実行可能解を求めるのは容易な場合が多い
整数計画問題
：実行可能解を求めることも一般にＮＰ完全
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、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
整数計画について
① 整数条件による高い表現能力
⇒ 整数計画問題は様々な組合せ最適化問題を含み、応用例も多い
巡回セールスマン問題、スケジューリング問題、．．．
⇒ 整数計画問題が解けることはインパクトが大きい
② ただ、整数計画問題は解くのが困難であった
 組合せ最適化問題のように、問題固有の構造を持たない
 組合せ最適化問題：実行可能解を求めるのは容易な場合が多い
整数計画問題
：実行可能解を求めることも一般にＮＰ完全
⇒ 組合せ最適化問題ごとに解法を開発するのが一般的であった
③ しかし、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
定式化の例１：数独
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
数独の定式化
方法１
xij : (i, j)マスに入る数字（つま
り1  xij  9）
方法２
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
定式化：方法１
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
定式化：方法１
注１
 xij  xi  x   x 
 


b  x  8b
 
xij  xi  1  b  x   8b 
b   b   1

b  , b   0 または1

注２
constraint programmingの分野では
all_different( xi1,, xin )
という表現をする
53
定式化：方法２
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
定式化の比較
 定式化１
 変数 3970（一般整数変数 2025、0-1変数 1944）
 制約 4860＋固定制約数
 計算時間 41.5秒（Intel Core2 Duo [email protected], メモリ2.0GB,
CPLEX 12.5）
※Kochの論文によると、CPLEX9.03では6時間経っても解けなかった
 定式化２
 変数 729（0-1変数 729）
 制約 346＋固定制約数
 計算時間 0.02秒
55
巡回セールスマン問題（TSP）
Traveling Salesman Problem
各都市をちょうど１回ずつ訪問し、
出発地点に戻ってくる閉路（巡回
路）の中で、ウェイト（距離、時間
等）の合計が最小となるものを求
めよ
兵庫県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 , jS
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)
出る辺は１本
i
j
入る辺は１本
部分巡回路
除去制約
58
部分巡回路除去制約
Dantzig-Fulkerson-Johnson (1954) 再掲
 xij  | S | 1 (S  V , 1 S , | S | 1)
i , jS
 変数 𝑂 𝑛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
変数
制約式
初期ＬＰ
ＩＰ最終解
目的関数値
時間（秒）
上界値
下界値
誤差
時間（秒）
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
おわりに
まとめ１
 整数計画問題：「解けない」から「解ける」へ
MIP：汎用的 → インパクトが大きい
 現在のMIPソルバー = LPベースの分枝カット法
⇒ 苦手な問題の種類もある（big-Mを含む問題など）
まとめ２
 問題の捉え方によって定式化が異なる
 変数の設定方法によって定式化が異なる
⇒「悪い定式化」になる可能性も
 しかし、まずは定式化をしてみる
 安心のため．．．
 強引な定式化でも、今のソルバーなら解けてしまうかも
67