制約とペナルティ

制約付き最適化問題

数理最適化において、制約とは解が満たさなければならない条件のことです。例えば、次の問題は制約付き最適化問題になります。

\[\begin{split} \begin{aligned} \text{Minimize} & \quad f(x) \\ \text{subject to} & \quad g(x) = 0 \end{aligned} \end{split}\]

ここで\(f\)と\(g\)は決定変数\(x\)の関数です。条件\(g(x) = 0\)は等式制約と呼ばれます。\(g(x) = 0\)を満たすすべての\(x\)の集合は実行可能集合と呼ばれます。制約は\(g(x) = 0\)のような等式制約に限らず、\(g(x) \leq 0\)のような不等式制約である場合もあります。例えば、次の問題も制約付き最適化問題です。

\[\begin{split} \begin{aligned} \text{Minimize} & \quad f(x) \\ \text{subject to} & \quad g(x) \leq 0 \end{aligned} \end{split}\]

jijmodelingでは、Constraintクラスを使用して等式制約と不等式制約の両方を記述できます。例えば、等式制約\(\sum_i x_i = 1\)は次のように表現できます。

import jijmodeling as jm

N = jm.Placeholder("N")
x = jm.BinaryVar("x", shape=(N,))
n = jm.Element('n', belong_to=(0, N))

jm.Constraint("onehot", jm.sum(n, x[n]) == 1)

\[\begin{split}\begin{array}{cccc} & \text{onehot} & \displaystyle \sum_{n = 0}^{N - 1} x_{n} = 1 & \\ \end{array}\end{split}\]

上記のコードでは、jm.Constraintの第1引数として文字列”onehot”を指定していることに注意してください。制約オブジェクトには名前と制約式があります。これらの名前は制約が満たされているかどうかを確認する際に使用されます。制約式は、==、<=、または>=の3つの比較演算子のいずれかを使用する論理式でなければなりません。また、次のようにして1つの問題に複数の制約を課すことができます。

x = jm.BinaryVar("x", shape=(4,))
problem = jm.Problem("constraint_sample")
problem += jm.Constraint("c01", x[0] + x[1] <= 1)
problem += jm.Constraint("c12", x[1] + x[2] <= 1)
problem += jm.Constraint("c23", x[2] + x[3] >= 0)

Tip

他の比較演算子（例えば>）や論理演算子はサポートされていません。

x = jm.BinaryVar("x", shape=(4,))
jm.Constraint("unsupported", (x[0] + x[1] <= 1) | (x[1] + x[2] <= 1))

forall制約

しばしば制約は変数によってインデックス付けされます。例えば、次のような問題は制約付き最適化問題です。

\[\begin{split} \begin{aligned} \text{Minimize} & \quad \sum_{i=0}^{N-1} \sum_{j=0}^{M-1} a_{ij} x_{ij} \\ \text{subject to} & \quad \sum_{j = 0}^{M - 1} x_{ij} = 1 \quad \forall i \in \left\{0, \ldots, N - 1\right\} \end{aligned} \end{split}\]

このような\(\forall i \in \left\{0, \ldots, N - 1\right\}\)を表現するために、Constraintオブジェクトにはforallオプションがあります。例えば、上記の問題は次のように表現できます。

N = jm.Placeholder("N")
M = jm.Placeholder("M")
a = jm.Placeholder("a", ndim=2)
x = jm.BinaryVar("x", shape=(N, M))
i = jm.Element('i', belong_to=(0, N))
j = jm.Element('j', belong_to=(0, M))

problem = jm.Problem("forall_sample")
problem += jm.sum([i, j], a[i, j] * x[i, j])
problem += jm.Constraint("onehot", jm.sum(j, x[i, j]) == 1, forall=i)

ペナルティとは

ペナルティ法とラグランジュ乗数法 は制約付き最適化問題を制約なし最適化問題に変換するための最も一般的な方法です。ここではペナルティ法について見ていきましょう。

\[\begin{split} \begin{aligned} \text{Minimize} & \quad f(x) \\ \text{subject to} & \quad g(x) = 0 \end{aligned} \end{split}\]

この問題は次のような制約なし最適化問題に変換されます。

\[ \text{Minimize} \quad f(x) + \alpha p(x), \]

この変換では\(\alpha\)（ペナルティ係数またはラグランジュ乗数）と\(p(x)\)（ペナルティ項）が重要な役割を果たします。通常、\(p(x)\)は\(p(x) = g(x)^2\)として定義されます。\(f(x) + \alpha p(x)\)の最小値が\(p(x) = 0\)を満たす場合、その\(x\)は元の制約付き最適化問題の最小値になります。ペナルティ\(p(x)\)が正の場合、ペナルティ係数\(\alpha\)の値を増やして上記の制約なし最適化問題を解くと元の最適化問題の解を得られる可能性が高まります。

一部のソルバーは制約なし最適化問題のみしか受け入れません。QUBOの「U」は「Unconstrained（無制約）」を意味しています。jijmodelingで制約付き最適化問題として定式化した数理モデルを入力形式がQUBOのソルバーに与えるためには、jijmodeling_transpilerまたはJijZeptを利用して制約なし最適化問題に変換する必要があります。しかし、実際の最適化問題においては\(p\)の定義も重要であり、そのため、jijmodelingではペナルティ項をカスタマイズする機能を提供しています。

制約をペナルティ項に変換する

jijmodelingは制約をペナルティ項に変換する機能を持ちません。ここではjijmodeling_transpilerが制約をペナルティ項に変換する方法を説明します。以下のような簡単な問題を考えてみましょう。

\[\begin{split} \begin{aligned} \text{Minimize} & \quad \sum_{i=0}^{N-1} a_i x_i \\ \text{subject to} & \quad \sum_{i = 0}^{N - 1} x_i = 1 \end{aligned} \end{split}\]

この問題はjijmodelingで次のように定式化できます。

import jijmodeling as jm

a = jm.Placeholder("a", ndim=1)
N = a.len_at(0, latex="N")

i = jm.Element('i', belong_to=(0, N))
x = jm.BinaryVar('x', shape=(N,))

problem = jm.Problem('translate_constraint_to_penalty')
problem += jm.sum(i, a[i]*x[i])
problem += jm.Constraint(
    'onehot',
    jm.sum(i, x[i]) == 1,
)
problem

\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{translate\_constraint\_to\_penalty} & & \\& & \min \quad \displaystyle \sum_{i = 0}^{N - 1} a_{i} \cdot x_{i} & \\\text{{s.t.}} & & & \\ & \text{onehot} & \displaystyle \sum_{i = 0}^{N - 1} x_{i} = 1 & \\\text{{where}} & & & \\& x & 1\text{-dim binary variable}\\\end{array}\end{split}\]

jijmodeling_transpilerでは、この制約付き最適化問題は次の制約なし最適化問題に変換されます。

\[ \text{Minimize} \quad \sum_{i=0}^{N-1} a_i x_i + \alpha \left(\sum_{i = 0}^{N - 1} x_i - 1\right)^2 \]

ここでは、\(a = [1, 2]\)とし、\(\alpha = 5\)とする場合を見てみましょう。

import jijmodeling_transpiler as jmt

instance_data = {
    "a": [1, 2],
}

compiled_model = jmt.core.compile_model(problem, instance_data)
pubo_builder = jmt.core.pubo.transpile_to_pubo(
    compiled_model,
    normalize=False  # 簡単のために正規化を無効化
)
qubo, constant = pubo_builder.get_qubo_dict(multipliers={ 'onehot': 5 })

結果は次のようになります。

qubo

{(0, 0): -4.0, (1, 1): -3.0, (0, 1): 10.0}

constant

5.0

なぜ、このようなquboとconstantが得られるのかというと、以下のような計算が行われたからです。

\[\begin{split} \begin{aligned} \sum_{i=0}^{N-1} a_i x_i + \alpha \left(\sum_{i = 0}^{N - 1} x_i - 1\right)^2 &= x_1 + 2 x_2 + 5 (x_1 + x_2 - 1)^2 \\ &= -4 x_1 - 3 x_2 + 10 x_1 x_2 + 5 \\ &= \begin{bmatrix} x_1 & x_2 \end{bmatrix} \begin{bmatrix} -4 & 10 \\ 0 & -3 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} + 5 \end{aligned} \end{split}\]

上記の計算は、バイナリ変数\(x_i\)が\(x_i^2 = x_i\)であることを使って式変形しています。

jijmodeling_transpilerの変換過程には2つのフェーズに分かれています。

• Problemオブジェクトとinstance_dataを使用してCompiledInstanceオブジェクトに変換する。

• CompiledInstanceオブジェクトに未定乗数を指定してQUBOに変換する。

Note

jijmodeling_transpilerでは、get_qubo_dictメソッドのmultipilers引数を使用することで各ペナルティ項の未定乗数を設定することができます。また、jijmodeling_transpilerでは拡張ラグランジュ法のような他の緩和方法もサポートされています。詳細についてはjijmodeling_transpilerのリファレンスを確認してください。

CustomPenaltyTerm

制約をペナルティ項に変換する作業はjijmodeling_transpilerやJijZeptの各種サービスの役割ですが、必要に応じてオリジナルのペナルティ項を設定したい場合も存在します。jijmodelingではCustomPenaltyTermを使用することでペナルティ項のカスタマイズが可能になります。ここでは、前述の例と同じペナルティ項をCustomPenaltyTermで定義する方法を説明します。

\[ \text{Minimize} \quad \sum_{i=0}^{N-1} a_i x_i + \alpha \left(\sum_{i = 0}^{N - 1} x_i - 1\right)^2 \]

この問題は次のようにCustomPenaltyTermを使用して表現できます。

import jijmodeling as jm

a = jm.Placeholder("a", ndim=1)
N = a.len_at(0, latex="N")

i = jm.Element('i', belong_to=(0, N))
x = jm.BinaryVar('x', shape=(N,))

problem = jm.Problem('penalty_sample')
problem += jm.sum(i, a[i]*x[i])
problem += jm.CustomPenaltyTerm(
    'onehot',
    (jm.sum(i, x[i]) - 1)**2,
)

この時点では未定乗数\(\alpha\)が設定されていないことに注意してください。この場合においても前述と同じようにjijmodeling_transpilerを利用することができます。

import jijmodeling_transpiler as jmt

instance_data = {
    "a": [1, 2],
}

compiled_model = jmt.core.compile_model(problem, instance_data)
pubo_builder = jmt.core.pubo.transpile_to_pubo(
    compiled_model,
    normalize=False  # 簡単のために正規化を無効化
)
qubo, constant = pubo_builder.get_qubo_dict(multipliers={ 'onehot': 5 })

当然、結果は同じになります。

qubo

{(0, 0): -4.0, (1, 1): -3.0, (0, 1): 10.0}

constant

5.0

また、CustomPenaltyTermにもforall引数があります。

import jijmodeling as jm

a = jm.Placeholder("a", ndim=2)
N = a.len_at(0, latex="N")
M = a.len_at(1, latex="M")

i = jm.Element('i', belong_to=(0, N))
j = jm.Element('j', belong_to=(0, M))
x = jm.BinaryVar('x', shape=(N, M))

problem = jm.Problem('forall_penalty_sample')
problem += jm.sum([i, j], a[i, j]*x[i, j])
problem += jm.CustomPenaltyTerm(
    'onehot',
    (jm.sum(j, x[i, j]) - 1)**2,
    forall=i
)
problem

\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{forall\_penalty\_sample} & & \\& & \min \quad \displaystyle \sum_{i = 0}^{N - 1} \sum_{j = 0}^{M - 1} a_{i, j} \cdot x_{i, j} & \\\text{{penalty terms}} & & & \\ & \text{onehot} & \displaystyle \left(\left(\sum_{j = 0}^{M - 1} x_{i, j} - 1\right)^{2}\right) & \forall i \in \left\{0,\ldots,N - 1\right\} \\\text{{where}} & & & \\& x & 2\text{-dim binary variable}\\\end{array}\end{split}\]