この記事について
整数計画で困ることのひとつが「非線形なものを扱いたい」という場面です。区分関数($\max, \min$)、非凸な非線形関数、2次の積項…。これらを線形のMIPで表現できれば、既存のソルバー(Gurobi、CPLEX、HiGHS等)の強力なエンジンに乗せられます。
この記事では、区分線形近似と2値変数の積の線形化という2つの代表的なテクニックを紹介します。
非凸な非線形関数の区分線形近似
非凸な非線形関数 $f(x)$ を最小化したい場合、そのままだと整数計画では扱えません。そこで 区分線形関数 $g(x)$ で近似します。
ブレークポイント $a_1, a_2, \ldots, a_m$ を選び、その点での値 $f(a_i)$ を結んだ折れ線で近似します。ブレークポイントが多いほど近似精度は上がりますが、変数も増えるので計算が重くなります。
線分上の点の表現
区分線形関数の任意の点は「いずれかの線分上にある」ので、その線分の2端点の重み付き平均で書けます。
(x, g(x)) = t_i (a_i, f(a_i)) + t_{i+1} (a_{i+1}, f(a_{i+1}))
t_i + t_{i+1} = 1, \quad t_i, t_{i+1} \geq 0
これを全体に拡張すると、ブレークポイント $a_1, \ldots, a_m$ の重みベクトル $(t_1, \ldots, t_m)$ で:
(x, g(x)) = \sum_{i=1}^{m} t_i (a_i, f(a_i)), \quad \sum_{i=1}^{m} t_i = 1, \quad t_i \geq 0
$t_i$ のうち高々2つの隣り合うものだけが正でなければなりません。これがないと、隣同士でない2点の重み付き平均になってしまい、関数 $f$ と値がずれます。
「隣り合う2つだけが正」を整数計画で表現
各線分 $i$ に対応する2値変数 $z_i \in {0, 1}$ を導入し、「点が $i$ 番目の線分上にある」ことを表します。
\begin{aligned}
t_1 &\leq z_1 \\
t_i &\leq z_{i-1} + z_i, \quad i = 2, \ldots, m-1 \\
t_m &\leq z_{m-1} \\
\sum_{i=1}^{m-1} z_i &= 1 \\
z_i &\in \{0, 1\}
\end{aligned}
- $z_i = 1$ なら点は $i$ 番目の線分上にある($t_i$ と $t_{i+1}$ のみが正になり得る)
- $\sum z_i = 1$ により、ちょうど1つの線分が選ばれる
これは SOS2(Special Ordered Set type 2) という制約タイプとして Gurobi・CPLEX に標準搭載されています。直接 SOS2 制約として渡せば、より効率的に処理してくれます。
Pythonで区分線形近似を実装
import numpy as np
import pulp
# 近似したい非線形関数(例:sin に少し小細工を加えた非凸関数)
def f(x):
return np.sin(x) + 0.1 * x
# ブレークポイント
xs = np.linspace(0, 10, 11) # 11点
fs = [f(x) for x in xs]
m = len(xs)
prob = pulp.LpProblem("piecewise", pulp.LpMinimize)
# 区分線形近似に使う変数
t = [pulp.LpVariable(f"t{i}", lowBound=0, upBound=1) for i in range(m)]
z = [pulp.LpVariable(f"z{i}", cat='Binary') for i in range(m-1)]
# x と g(x) を t で表現
x_var = pulp.lpSum(xs[i] * t[i] for i in range(m))
g_x = pulp.lpSum(fs[i] * t[i] for i in range(m))
# 目的関数:g(x) を最小化
prob += g_x
# SOS2 制約(隣り合う2つの t のみが正)
prob += t[0] <= z[0]
for i in range(1, m-1):
prob += t[i] <= z[i-1] + z[i]
prob += t[m-1] <= z[m-2]
prob += pulp.lpSum(t) == 1
prob += pulp.lpSum(z) == 1
prob.solve(pulp.PULP_CBC_CMD(msg=0))
x_opt = pulp.value(x_var)
g_opt = pulp.value(g_x)
print(f"近似解 x = {x_opt:.4f}")
print(f"近似 g(x) = {g_opt:.4f}")
print(f"真の f(x) = {f(x_opt):.4f}")
実行結果例:
近似解 x = 5.0000
近似 g(x) = -0.4589
真の f(x) = -0.4589
これで非凸関数の大域最適解に近い点を見つけられます。
2値変数の積を線形化する
2次計画問題(QP)や QUBO(Quadratic Unconstrained Binary Optimization)を MIP に変換するときに重要な技法です。
問題
\text{minimize} \quad \sum_{i=1}^{n} \sum_{j=1}^{n} q_{ij} x_i x_j
\text{subject to} \quad x_i \in \{0, 1\}
$x_i x_j$ の積項が非線形です。これを線形化したいです。
古典的な線形化
新しい2値変数 $y_{ij} \in {0, 1}$ を導入して $y_{ij} = x_i x_j$ を表します。$x_i, x_j \in {0,1}$ の組み合わせは4通り:
| $x_i$ | $x_j$ | $x_i x_j$ |
|---|---|---|
| 0 | 0 | 0 |
| 1 | 0 | 0 |
| 0 | 1 | 0 |
| 1 | 1 | 1 |
これを満たすには3本の線形制約で十分です。
\begin{aligned}
y_{ij} &\leq x_i \\
y_{ij} &\leq x_j \\
y_{ij} &\geq x_i + x_j - 1 \\
y_{ij} &\in \{0, 1\}
\end{aligned}
確認
| $x_i$ | $x_j$ | 上限 $\min(x_i, x_j)$ | 下限 $x_i + x_j - 1$ | $y_{ij}$ |
|---|---|---|---|---|
| 0 | 0 | 0 | -1 | 0 |
| 1 | 0 | 0 | 0 | 0 |
| 0 | 1 | 0 | 0 | 0 |
| 1 | 1 | 1 | 1 | 1 |
線形化後の MIP
\begin{aligned}
\text{minimize} \quad & \sum_{i=1}^{n} \sum_{j=1}^{n} q_{ij} y_{ij} \\
\text{subject to} \quad & y_{ij} \leq x_i \\
& y_{ij} \leq x_j \\
& y_{ij} \geq x_i + x_j - 1 \\
& x_i \in \{0, 1\} \\
& y_{ij} \in \{0, 1\}
\end{aligned}
これで標準的なMIPソルバーで解けます。
Pythonでの実装:QUBO → MIP変換
import numpy as np
import pulp
# QUBO 行列の例(対称・非対角に交互作用あり)
np.random.seed(42)
n = 5
Q = np.random.randint(-5, 5, (n, n))
Q = (Q + Q.T) / 2 # 対称化
prob = pulp.LpProblem("qubo_to_mip", pulp.LpMinimize)
x = [pulp.LpVariable(f"x{i}", cat='Binary') for i in range(n)]
y = {(i, j): pulp.LpVariable(f"y{i}_{j}", cat='Binary')
for i in range(n) for j in range(n)}
# 線形化された目的関数
prob += pulp.lpSum(Q[i][j] * y[(i, j)] for i in range(n) for j in range(n))
# 線形化制約
for i in range(n):
for j in range(n):
prob += y[(i, j)] <= x[i]
prob += y[(i, j)] <= x[j]
prob += y[(i, j)] >= x[i] + x[j] - 1
prob.solve(pulp.PULP_CBC_CMD(msg=0))
x_opt = [int(pulp.value(xi)) for xi in x]
print(f"x* = {x_opt}")
print(f"目的関数値 = {pulp.value(prob.objective):.2f}")
# 検証:直接 xQx を計算
x_np = np.array(x_opt)
direct = x_np @ Q @ x_np
print(f"直接計算 xQx = {direct:.2f}")
実行結果例:
x* = [0, 1, 0, 1, 1]
目的関数値 = -12.00
直接計算 xQx = -12.00
線形化した目的関数値と直接計算した $\mathbf{x}^\top Q \mathbf{x}$ が一致しています。これで非線形の2次問題を線形のMIPで解けました。
まとめ
- 区分線形近似:非凸な非線形関数 $f(x)$ をブレークポイントで折れ線近似
- 「隣り合う2つの重みのみが正」を SOS2 制約 で表現
- 2値変数の積の線形化:$y_{ij} = x_i x_j$ を3本の線形制約で表現
- これによりQUBO(2次の0-1最適化)が MIP として解ける
