混合整数計画問題（MIP）定式化チートシート（PuLP）

PuLPソルバーで解きたいけど, 式変形がよくわからないという人向けのチートシート.
導出する機会があったため, メモとして残しておく.
(指示関数の導出に苦労したので記事を書く. 他はおまけ. )

Indicator function（指示関数）

条件満たすとき1, そうでないとき0を表したい.

$$
\mathbb{I}[x\le a]=\begin{cases}
1&, x\le a\\
0&, x > a
\end{cases}
$$

前提として, 次の表で登場する変数$x$は定数$m,M$を用いて以下を満たす.
$$ m \le x \le M $$
条件をみたすかどうかを表すバイナリ変数を$d \in \{0, 1\}$とする.
等号の無い不等号$<,>$は, PuLPで利用できないため, 小さな正値$0<\varepsilon\ll 1$を使って表現する.

表現したいこと	定式化	確認
$$x \le a \to d=1\\ \Updownarrow \text{ (equiv.)} \\ d=0 \to a < x $$	$x - a \ge (m-a)d + \varepsilon (1-d)$	$ \begin{split} & d=0\to x \ge a+\varepsilon \to x > a, \\ &d=1 \to x \ge m \text{ (always true)}\end{split}$
$$a < x \to d=0\\ \Updownarrow \text{ (equiv.)} \\ d=1 \to x \le a $$	$x - a \le (M-a)(1-d)$	$\begin{split} & d=1\to x \le a, \\ &d=0\to x \le M\text{ (always true)}\end{split}$
$$\mathbb{I}[ x \le a ]$$	$(m-a)d + \varepsilon (1-d) \le x - a \le (M-a)(1-d)$	上2つをまとめただけ
$$\mathbb{I}[ y \ge b ]$$	$(b-M’)d + \varepsilon (1-d) \le b - y \le (b-m’)(1-d)$	$$\text{let } y=-x, b=-a, m’=-M, M’=-m, \\\ \text{then } m’ \le y \le M’$$

$A\Leftrightarrow B$は同値を表しているので, $A$と$B$は同じ定式化になる.
$A$と$B$は対偶を取っただけですが, 証明のしやすさが違います.

指示関数を使って論理式を作り, AND($d_1=1, d_2=1$)やOR($d_1 + d_2\ge 1$)などを表現することもできる.

その他の式変形

以下では, $x,y,z,w$を変数, $a,b,n,M$を定数とする.

表現したいこと	定式化	補足
$$\min_x |x| $$	$$ \begin{split} \min_{x,y} ~&~ y \\ \mathrm{s.t.}~&~x\le y, -x\le y \end{split} $$
$$y=|x|$$	$$y^+-y^-=x, \\ y^++y^-=y, \\ y^+, y^-\ge 0,\\ z\in\{0,1\}, \\ y^+\le Mz,\\y^-\le M(1-z) $$	$M$は$\max|x|$以上の定数. 変数$y$は, $x$が正のとき$y^+$で表され, 負のとき$y^-$で表現される.
$$ \min_{x,y} \max\{x, y\} $$	$$\begin{split}\min_{x,y,z}~&~ z\\ \mathrm{s.t.}~&~x\le z, y\le z \end{split}$$	$y=0$のケースはよく使う
$$x=a\text{ or }x=b$$	$$y\in\{0,1\}, \\ x=ya+(1-y)b$$
$$ x,y\in\mathbb{Z}, \\ 1\le x \le n, 1\le y \le n\\ x\ne y $$	$$ z,w\in \{0,1\}^n,\\ \sum_i z_i = 1, \sum_i w_i = 1, \\ z_i+w_i\le 1, \forall i $$	整数値をonehotベクトルで表現して, 2つの整数が異なることを和が1以下として表す. 数独の定式化もこれ.
$$\min_{x} \frac{x^\top a}{x^\top b} $$	$$\begin{split}\min_{y} ~&~y^\top a, \\ \mathrm{s.t.}~&~y^\top b=1\end{split}$$

更に複雑な変形が必要なら, "整数計画 定式化 テクニック"や"mip formulation technique"などで調べるとよい.
追記: こんなページも見つけた→定式化技法集

簡単にできる高速化

スレッド数指定

problem.solve(pulp.PULP_CBC_CMD(msg=0, threads=4))

• threadsで並列数を指定できる
• msg=0でpulpによる標準出力を表示させないようにすることができる

連続変数の代わりに整数変数にしてみる

x = LpVariable("x", 0, 1, pulp.LpContinuous)
# 上記の代わりに
PRECISION = 1000
x = LpVariable("x", 0, PRECISION, pulp.LpInteger) * (1.0 / PRECISION)

• 細かい連続値の値が不要なときに使える
• 速くなるかは問題依存な気がする

おまけ（検算）

PuLPで指示関数の挙動を確認してみる

import pulp

# 小さな値
eps = 1e-6
# xの範囲
m = 1
M = 10
# 変数
x = pulp.LpVariable("x", m, M, cat="Integer")
d = pulp.LpVariable("d", cat="Binary")
# 最大化問題
problem = pulp.LpProblem("problem", pulp.LpMaximize)
problem += x
# d = I[x <= a] を確認する
a = 3
problem += x - a >= (m-a)*d + eps*(1-d)
problem += x - a <= (M-a)*(1-d)
# 以下のどれかをonにして, 挙動を確認できる
# problem += x == 5 # -> d = 0
# problem += x == 2 # -> d = 1
# problem += d == 1 # -> x = a
# problem += d == 0 # -> x = M
# 解く
problem.solve()
# ステータス
print(pulp.LpStatus[problem.status])
# 解
print(pulp.value(x), pulp.value(d))

コメントアウトしてある箇所を外して挙動を確認できます.