混合整数計画問題（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以下として表す. [数独の定式化](https://qiita.com/birdwatcher/items/d5c2a88c733ca71ae91e#%E5%AE%9A%E5%BC%8F%E5%8C%96)もこれ. |
| $$\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"などで調べるとよい. 
追記: こんなページも見つけた→[定式化技法集](https://www.msi.co.jp/solution/nuopt/docs/techniques/index.html) 

# 簡単にできる高速化
**スレッド数指定**
```python
problem.solve(pulp.PULP_CBC_CMD(msg=0, threads=4))
```
- `threads`で並列数を指定できる
- `msg=0`でpulpによる標準出力を表示させないようにすることができる

**連続変数の代わりに整数変数にしてみる**
```python
x = LpVariable("x", 0, 1, pulp.LpContinuous)
# 上記の代わりに
PRECISION = 1000
x = LpVariable("x", 0, PRECISION, pulp.LpInteger) * (1.0 / PRECISION)
```
- 細かい連続値の値が不要なときに使える
- 速くなるかは問題依存な気がする

# おまけ（検算）
PuLPで指示関数の挙動を確認してみる
```python
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))
```

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