PuLP での論理和と論理積の線形化

前回の記事では補助変数を利用することで、論理包含のある一見すると非線形な制約を線形制約へ変形し、PuLP で整数制約問題を解くことができました。

現実の問題を扱う際には、「ある特定の状況が望ましくないので、この条件を満たす場合にのみペナルティを足そう」というモデリングを行いたくなることが容易に想像できます。

そこで、論理積の線形表現を PuLP で実現する抽象を導出してみました。論理和もついでにやっておきました。

二項論理和 (OR) の線形化

論理和の方が簡単なので、先にやります。考え方は真理値表をそのまま実現することです。

$y = a \lor b$ を実現するために、次のように制約を置きます。

$y \geq a$ と $y \geq b$ により $a,b$ どちらかが $1$ であれば $y$ は $1$ になります。 $y \leq a + b$ は $a,b$ がともに $0$ のときに、 $y$ を $0$ にするための制約です。

$y \leq a + b$ の右辺の最大値は $2$ ですが、 $y$ は $0$ か $1$ なので $1$ で止まります。

このように制約付けることで、上真理値表のとおり唯一の $y$ が決まり $y = a \lor b$ となります。

二項論理積 (AND) の線形化

OR より難しいですが、考え方は同じです。

$z \leq a$ と $z \leq b$ で偽値があるときに $0$ となるように上から抑え、そうでないときは $z \geq a + b - 1$ によって $1$ にします。

これで $z = a \land b$ となります。

定数 $p$ を決めて、目的関数に $p \times z$ を足してやれば、「 $a$ と $b$ が両方 $1$ ならペナルティ $p$ 」を表現できます。

一般化

n 項の OR

n 項の AND

PuLP 実装サンプル

適当な最小化問題を作成してみました。 1 をとる変数を増やすとコストが減りますが、全部 1 だと大きなペナルティを食らう目的関数になっています。

import pulp as pl

prob = pl.LpProblem("DEMO", pl.LpMinimize)

# 既存の決定変数
a = pl.LpVariable("a", lowBound=0, upBound=1, cat="Binary")
b = pl.LpVariable("b", lowBound=0, upBound=1, cat="Binary")
c = pl.LpVariable("c", lowBound=0, upBound=1, cat="Binary")


# OR に条件付けするユーティリティ
def choreograph_or(prob, dest_var, *in_vars):
    for v in in_vars:
        prob += dest_var >= v
    prob += dest_var <= pl.lpSum(in_vars)

# AND に条件付けするユーティリティ
def choreograph_and(prob, dest_var, *in_vars):
    n = len(in_vars)
    for v in in_vars:
        prob += dest_var <= v
    prob += dest_var >= pl.lpSum(in_vars) - (n - 1)

# OR(a,b,c) を表す補助変数
y = pl.LpVariable("y", lowBound=0, upBound=1, cat="Binary")
choreograph_or(prob, y, a, b, c)

# AND(a,b,c) を表す補助変数
z = pl.LpVariable("z", lowBound=0, upBound=1, cat="Binary")
choreograph_and(prob, z, a, b, c)

# 目的関数（真なる変数を増やすとコストが減るが、全部真だと大きなペナルティを食らう）
prob += 5 * z - 2 * y - a - b - c


prob.solve()

print(
    f"{pl.value(a)=}",
    f"{pl.value(b)=}",
    f"{pl.value(c)=}",
    f"{pl.value(y)=}",
    f"{pl.value(z)=}",
    f"{pl.LpStatus[prob.status]=}",
    f"{pl.value(prob.objective)=}",
    sep="\n",
)

出力

pl.value(a)=0.0
pl.value(b)=1.0
pl.value(c)=1.0
pl.value(y)=1.0
pl.value(z)=0.0
pl.LpStatus[prob.status]='Optimal'
pl.value(prob.objective)=-4.0

AND が機能し a, b, c すべてを 1 にしないよう最適化されたことが分かります。

SCIP の場合

ところで SCIP では、AND/OR/XOR 制約などをそのまま表現でき、これらの制約に対応する効率的な最適化メソッドが用意されています¹。

model.addConsAnd([a, b, c], z)

しかし、目的関数を非線形にすることは SCIP でもできない²ので、補助変数にコンディションを表現させて係数をかけるというやり方は SCIP においても使います。