制約条件とペナルティ関数#

Amplify SDK では、任意の変数と多項式の次数を持つ制約条件を含むモデルを作成できます。しかし、組み合わせ最適化ソルバーごとに、制約条件として扱える変数の種類や次数は異なり、特に QUBO ソルバーでは制約条件そのものを受け取れないものもあります。

Amplify SDK は、目的関数と同様に、制約条件の変数と多項式の次数を変換することで、組み合わせ最適化ソルバーに対応した中間モデルを構築します。もし、このような変換で対応ができない場合には、制約条件に対応する ペナルティ関数 を生成し、自動で目的関数にペナルティ関数を足し込むことでソルバーを呼び出します。これらの処理により、制約条件をそのままでは受け取れない組み合わせ最適化ソルバーに対しても、制約条件を含むモデルの求解を可能にします。

中間モデルにおける制約条件#

それぞれのソルバークライアントには、目的関数に加えて 等式制約 と 不等式制約 それぞれについて、扱える変数の種類と次数が定義されています。Amplify SDK は、これらの情報をもとに、中間モデルにおける制約条件の変数と多項式の次数を決定して変換を試みます。

整数変数に対する等式制約を含む入力モデルを中間モデルに変換する例をいくつか紹介します。

from amplify import (
    VariableGenerator,
    Model,
    FixstarsClient,
    AcceptableDegrees,
    equal_to,
)

gen = VariableGenerator()
n = gen.scalar("Integer", bounds=(-10, 10))  # 整数変数を発行
c = equal_to(n, 1)  # n = 1 の等式制約を作成

model = Model(c)

例えば一次のバイナリ変数からなる等式制約をそのまま扱えるソルバークライアント向けには、次のように入力モデルから中間モデルへの変換が行われます。

bqbl = AcceptableDegrees(
    objective={"Binary": "Quadratic"},
    equality_constraints={"Binary": "Linear"}
)
im, mapping = model.to_intermediate_model(bqbl)

>>> print(im)
minimize:
  0
subject to:
  q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10 == 1 (weight: 1)
>>> print(mapping[n])
q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10

目的関数における変数変換と同様に、中間モデルに含まれる制約条件がバイナリ変数のみで構成され、中間モデルにおいて整数変数とバイナリ変数の変数変換が行われていることがわかります。

一方で、制約を全く扱えないソルバーに対しても、Amplify SDK は制約条件を含むモデルの求解が可能です。次の例は、QUBO ソルバーに対して、整数変数を含む入力モデルを中間モデルに変換する例です。

bq = AcceptableDegrees(objective={"Binary": "Quadratic"})
im, mapping = model.to_intermediate_model(bq)

>>> print(im)
minimize:
  0
subject to:
  q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10 == 1 (weight: 1)
>>> print(mapping[n])
q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10

一見すると前の例と全く違いが無いように見えますが、実は次節で説明するペナルティ関数と呼ばれる多項式が生成され、さらに変数変換されています。ペナルティ関数は目的関数の一部として利用されるため、目的関数と同様の変数変換と次数下げが行われます。これは次のようにして確認できます。

>>> # 入力モデルの制約条件のペナルティ関数
>>> print(model.constraints[0].penalty)
n_0^2 - 2 n_0 + 1

>>> # 中間モデルの制約条件のペナルティ関数
>>> print(im.constraints[0].penalty)  
4 q_0 q_1 + 8 q_0 q_2 + 16 q_0 q_3 + 10 q_0 q_4 + 16 q_1 q_2 + 32 q_1 q_3
  + 20 q_1 q_4 + 64 q_2 q_3 + 40 q_2 q_4 + 80 q_3 q_4 - 21 q_0 - 40 q_1
  - 72 q_2 - 112 q_3 - 85 q_4 + 121

ペナルティ法#

Amplify SDK は、ソルバーが扱える制約条件の種類や次数に応じて、変数変換や次数下げによって制約条件の条件式の変換を試みます。しかし変換ができない場合には、制約条件の条件式の代わりにペナルティ関数と呼ばれる多項式を目的関数に付与します。ペナルティ関数は制約条件と等価と見なせるように生成されるため、制約条件が間接的に実現されることになります。中間モデルを構築する際には、必要に応じて目的関数に加えてペナルティ関数の変数変換と次数下げも行われます。

ペナルティ法

変数 \(x_1, x_1, \ldots, x_n\) についてのある制約条件 \(c\) に対し、実数値をとる関数 \(p\) が

\[\begin{split} p(x_1, x_2, \ldots, x_n) \begin{cases} = 0 \quad & \text{if } c \text{ is satisfied} \\ > 0 \quad & \text{otherwise} \end{cases} \end{split}\]

をみたすとき、\(p\) を ペナルティ関数 といいます。

組合せ最適化最適化問題

\[\begin{split} \begin{align*} \text{minimize} & \quad f(x) \\ \text{subject to} & \quad g_1(x) \leq c_1, \, g_2(x) \leq c_2, \, \ldots, \, g_m(x) \leq c_m \end{align*} \end{split}\]

に対し、\(m\) 個の制約条件のペナルティ関数 \(p_1, p_2, \ldots, p_m\) をそれぞれ計算し、制約条件を持たない組合せ最適化問題

\[ \text{minimize} \quad f(x) + k_1 p_1(x) + k_2 p_2(x) + \cdots + k_m p_m(x) \]

に変換する方法を ペナルティ法 といいます。ここで、\(k_1, k_2, \ldots, k_m\) は十分に大きなハイパーパラメータです。

Constraint クラスのインスタンスに対して、そのペナルティ関数は、penalty プロパティを用いて取得できます。

例として、equal_to() ヘルパー関数を用いて作成された制約条件オブジェクトのペナルティ関数は、次のようにして確認できます。

from amplify import VariableGenerator, equal_to

gen = VariableGenerator()
q = gen.array("Binary", 6)
c = equal_to(q[0] + q[1] + q[2], 1)

>>> print(c.penalty)
2 q_0 q_1 + 2 q_0 q_2 + 2 q_1 q_2 - q_0 - q_1 - q_2 + 1

c.penalty の値は、c がみたされるとき、つまり q[0] + q[1] + q[2] == 1 となるときに 0 をとり、それ以外の場合は 0 よりも大きな値をとるようになっています。つまり、前節で定義されたペナルティ関数の要件をみたしています。

ペナルティ関数の重み#

ペナルティ関数は目的関数に追加されることで、制約を満たしていないときにのみ目的関数の値が大きくなる、まさにペナルティとして扱われます。ソルバーは目的関数とペナルティの値の合計を最小化するように解を探索するため、目的関数もペナルティも小さな解が得られることが期待されます。

ペナルティ法が上手く働くための重要な点として、制約条件の重み (上記の \(k_1\), \(k_2\), \(\ldots\), \(k_m\)) を適切に設定する必要があることが挙げられます。ソルバーは目的関数とペナルティの値の合計を最小化するため、制約条件を破ったときにかかるペナルティの値が小さい場合、制約条件を満たす中で目的関数が最小の解よりも、制約条件を満たさないが目的関数がとても小さい解の方が良い解だと判定される可能性があるためです。言い換えると、ソルバーに目的関数の最小化よりもペナルティ値の最小化を優先させる必要があるということになります。

ペナルティ関数の重みは、制約条件オブジェクト Constraint が持つ weight プロパティを用いて取得・設定できます。weight の初期値は 1 に設定されており、目的関数に足されるペナルティ関数は Amplify SDK が自動で計算する値に weight を乗じたものとなります。

gen = VariableGenerator()
q = gen.array("Binary", shape=(2, 3))
c_list = equal_to(q, 1, axis=1)

>>> print(c_list)
[q_{0,0} + q_{0,1} + q_{0,2} == 1 (weight: 1),
 q_{1,0} + q_{1,1} + q_{1,2} == 1 (weight: 1)]
>>> c_list[0].weight
1.0

次節で述べるように、通常は Amplify SDK が自動でペナルティ関数をペナルティとして最低 1 の値がかかるように規格化して設定します。つまり初期状態では 1 の値がペナルティとして課せられることになります。weight をどのくらいの値にするのが適切かは問題によって異なりますが、一つの目安としては、制約条件に関係する目的関数の値よりも大きく設定しておくと実行可能解が出やすくなります。そのため、その値が 1 からかけ離れている場合、ペナルティ関数の重みを適切に設定する必要があります。

ペナルティ関数の重みは直接 weight プロパティに代入して設定できますが、制約条件オブジェクト Constraint や制約条件リスト ConstraintList に対して数値をかけることでも設定できます。後者は複数の制約条件オブジェクトに一括で同じ値の重みを設定したいときや、モデルを構築した後で重みを設定したいときに便利です。

>>> c_list *= 2.0
>>> print(c_list)
[q_{0,0} + q_{0,1} + q_{0,2} == 1 (weight: 2),
 q_{1,0} + q_{1,1} + q_{1,2} == 1 (weight: 2)]

具体的な問題におけるペナルティ関数の重みの設定については下記を参照してください。

ヒント

たとえば次のような手順でペナルティ関数の重みを調整すると良いでしょう。

ある制約に着目したときにその制約が満たされていないときの 目的関数の最大の利得 を見積もります
- 例えば巡回セールスマン問題であれば、都市を 訪問しない ことによる最大の利得は、都市間の辺の長さの最大値と見積もれます
見積もった値よりも大きな値を対象の制約のペナルティ関数の重みとして設定します
実行可能解が得られなかった場合は重みを大きくしてソルバーを再度実行します
- 例えば 2 倍ずつ大きくしていくと良いでしょう

解の精度を上げたい場合には上記の後に以下の手順を行います。

重みを少しずつ小さくしてソルバーを複数回実行します
- 例えば 0.9 倍ずつ小さくしていくと良いでしょう
得られた全ての解のうち目的関数の値が最小のものを最終的な解とします

ペナルティ関数の自動生成#

equal_to() や less_equal() などのヘルパー関数を用いて制約条件オブジェクト Constraintを作成すると、Amplify は必要に応じて自動で最適なペナルティ関数を生成します。

Amplify SDK が制約条件のペナルティ関数を生成するとき、まず、制約条件式の左辺が取りうる値の上限と下限を大まかに見積もります。この見積もりは、たとえばバイナリ変数多項式の場合は、その多項式に含まれる負の係数の総和と正の係数の総和を計算することで得ることができます。

その後、等式・不等式制約の種類や指定されたアルゴリズムに応じて、以下のようにペナルティ関数を生成します。

等式制約#

等式制約 \(f(x) = c\) に対して、以下のようにペナルティ関数 \(p\) が生成されます。

(i) 制約条件式の左辺 \(f\) の取りうる値の範囲の下限が \(c\) と一致する場合

ペナルティ関数 \(p\) を \(f - c\) とします。

例えばバイナリ変数の積 q[0] * q[1] = 0 に対しては、\(f\) の下限と右辺 \(c\) が一致しているので、以下のようなペナルティ関数が生成されます。

>>> c = equal_to(q[0] * q[1], 0)
>>> print(c.penalty)
q_0 q_1

(ii) 制約条件式の左辺 \(f\) の取りうる値の範囲の上限が \(c\) と一致する場合

ペナルティ関数 \(p\) を \(c - f\) とします。

例えばバイナリ変数の積 q[0] * q[1] = 1 に対しては、\(f\) の上限と右辺 \(c\) が一致しているので、以下のようなペナルティ関数が生成されます。

>>> c = equal_to(q[0] * q[1], 1)
>>> print(c.penalty)
- q_0 q_1 + 1

(iii) 上記以外の場合

ペナルティ関数 \(p\) を \((f - c)^2\) とします。

例えばバイナリ変数の和 q[0] + q[1] + q[2] = 2 に対しては、以下のようなペナルティ関数が生成されます。

>>> c = equal_to(q[0] + q[1] + q[2], 2)
>>> print(c.penalty)
2 q_0 q_1 + 2 q_0 q_2 + 2 q_1 q_2 - 3 q_0 - 3 q_1 - 3 q_2 + 4

注釈

上記は、上限値と下限値が一致している不等式制約など等式制約とみなせる場合にも適用されます。

不等式制約#

不等式制約のペナルティ関数を生成する場合、まず、制約条件式の左辺が取りうる値の上界と下界の値を用いて制約条件を \(a \leq f \leq b\) という形に書き直します。この形の不等式制約に対してペナルティ関数を生成するアルゴリズムとして、以下のアルゴリズムが提供されています。アルゴリズムの指定は、不等式制約を生成するヘルパー関数 less_equal()、greater_equal()、clamp() の penalty_formulation キーワード引数を用いて行います。デフォルトは Default です。

Default (Default)

制約条件式の変数と係数がすべて整数値の場合は IntegerVariable アルゴリズムを使用し、そうでない場合は RealVariable アルゴリズムを使用します。

IntegerVariable

整数値をとる補助変数を用いてペナルティを生成します。不等式制約 \(a \leq f \leq b\) に対して \(a\) 以上 \(b\) 以下の値を取る整数変数 \(n\) を発行し、等式制約 \(f - n = 0\) のペナルティ関数を不等式制約 \(a \leq f \leq b\) のペナルティ関数とします。

ただし、例外的に、\(b - a = 1\) となっている場合は、整数変数を発行せずに \((f - a)(f - b)\) がペナルティ関数となります。

注意

整数値でない制約条件に IntegerVariable を指定した場合、正確な定式化にはならないので注意してください。

例えばバイナリ変数の和 q[0] + q[1] + q[2] <= 2 に対しては、以下のようなペナルティ関数が生成されます。

>>> c = less_equal(q[0] + q[1] + q[2], 2, penalty_formulation="IntegerVariable")
>>> print(c.penalty)  
2 q_0 q_1 + 2 q_0 q_2 - 2 q_0 n_0 + 2 q_1 q_2 - 2 q_1 n_0 - 2 q_2 n_0
  + n_0^2 + q_0 + q_1 + q_2
>>> print(gen.variables[3])
{name: n_0, id: 3, type: Integer, lower_bound: -0, upper_bound: 2}

ここで、n_0 はペナルティ関数の生成と共に発行された整数の補助変数です。

注釈

QUBO ソルバーなど整数変数を扱えないソルバーでは、中間モデルの構築時に変数変換が行われます。詳細は「変数変換と次数下げ」を参照してください。

この変数変換では、整数変数の取り得る値の範囲が大きいほど必要な補助変数の数は多くなります。そのため、制約条件の構築時点で両辺を約分するなどして、不等式制約の取り得る範囲を小さくしておくと変数変換が効率化されることがあります。

バイナリ変数の和 q[0] + q[1] + q[2] <= 1 に対しては、制約条件の左辺 \(f\) の下限 \(0\) と右辺 \(c=1\) の差が 1 なので、整数変数は発行されず以下のようなペナルティ関数が生成されます。

>>> c = less_equal(q[0] + q[1] + q[2], 1, penalty_formulation="IntegerVariable")
>>> print(c.penalty)
2 q_0 q_1 + 2 q_0 q_2 + 2 q_1 q_2

RealVariable

実数値をとる補助変数を用いてペナルティを生成します。不等式制約 \(a \leq f \leq b\) に対して \(a\) 以上 \(b\) 以下の値を取る実数変数 \(x\) を発行し、等式制約 \(f - x = 0\) のペナルティ関数を不等式制約 \(a \leq f \leq b\) のペナルティ関数とします。

例えばバイナリ変数の和 0.1 * q[0] + 0.2 * q[1] + 0.4 * q[2] <= 0.5 に対しては、以下のようなペナルティ関数が生成されます。

>>> c = less_equal(0.1 * q[0] + 0.2 * q[1] + 0.4 * q[2], 0.5, penalty_formulation="RealVariable")
>>> print(c.penalty)
0.04 q_0 q_1 + 0.08 q_0 q_2 - 0.2 q_0 x_0 + 0.16 q_1 q_2 - 0.4 q_1 x_0 - 0.8 q_2 x_0 + x_0^2 + 0.01 q_0 + 0.04 q_1 + 0.16 q_2
>>> print(gen.variables[3])
{name: x_0, id: 3, type: Real, lower_bound: 0, upper_bound: 0.5}

ここで、x_0 はペナルティ関数の生成と共に発行された実数の補助変数です。

注意

現在の Amplify SDK では実数変数からバイナリ変数への変数変換はサポートされていません。そのため、QUBO ソルバーなど整数変数を扱えないソルバーでは、中間モデルの構築に失敗する可能性があります。

対処としては以下の2つが考えられます。

制約条件の両辺を定数倍し両辺を整数値にして IntegerVariable を使用する
以下で述べる緩和法 (Relaxation) による近似を試みる

Relaxation

制約条件の左辺 \(f\) が二次以上で、制約の範囲と \(f\) の下限あるいは上限が一致している場合は LinearRelaxation を使用し、そうでない場合は QuadraticRelaxation を使用します。

注意

緩和法 (Relaxation, LinearRelaxation, QuadraticRelaxation) によって作られるペナルティ関数は正確な定式化ではないことに注意してください。「ペナルティ法」で定義されたペナルティ関数の要件をみたさないためです。

LinearRelaxation が適用される場合は \(f\) を \(a\) あるいは \(b\) に近づける方向に、QuadraticRelaxation が適用される場合は \(f\) を \((a + b) / 2\) に近づける方向にペナルティがかかります。しかし、制約を満たす場合にかかるペナルティ値が一定でないため、ペナルティの重みを十分に大きくしても最適解が見つかるとは限りません。そのため、より良い解を見つけるためには、ペナルティ関数の重みを変更しながら何度もソルバーを実行し、実行可能解の中から最良解を選択するという方針が考えられます。

LinearRelaxation

ペナルティ法ではなくラグランジュ緩和を用います。
不等式制約 \(a \leq f \leq b\) に対して、\(f\) の取りうる値の範囲の下限が \(a\) と一致する場合、\(f - a\) を規格化したものを生成し、\(f\) の取りうる値の範囲の上限が \(b\) と一致する場合は、\(b - f\) を規格化したものを生成します。どちらでもない場合は QuadraticRelaxation が使用されます。

例として、バイナリ変数の和 q[0] + q[1] + q[2] <= 2 に対しては、以下のようなペナルティ関数が生成されます。

>>> c = less_equal(q[0] + q[1] + q[2], 2, penalty_formulation="LinearRelaxation")
>>> print(c.penalty)
0.5 q_0 + 0.5 q_1 + 0.5 q_2

QuadraticRelaxation

ペナルティ法ではなくラグランジュ緩和を用います。不等式制約 \(a \leq f \leq b\) に対して、\((f - \left(a + b\right)/2)^2\) を規格化したものを生成します。

例として、バイナリ変数の和 q[0] + q[1] + q[2] <= 2 に対しては、以下のようなペナルティ関数が生成されます。

>>> c = less_equal(q[0] + q[1] + q[2], 2, penalty_formulation="QuadraticRelaxation")
>>> print(c.penalty)
2 q_0 q_1 + 2 q_0 q_2 + 2 q_1 q_2 - q_0 - q_1 - q_2 + 1

ペナルティ関数の指定#

ヘルパー関数を用いて作成された制約条件に対しては、Amplify が自動でペナルティ関数を指定します。一方で、ユーザが自身で定義したペナルティ関数を設定したい場合には Constraint のコンストラクタを呼び出します。

次のように、制約条件の式とペナルティ関数を指定して制約条件オブジェクトを作成することができます。

制約条件式が \(q_0 + q_1 = 2\) でペナルティが \(- q_0 - q_1\) の不等式制約を作成する:
```
>>> c = Constraint(q[0] + q[1], eq=2, penalty=-q[0] - q[1])
```
制約条件式が \(q_0 + q_1 \leq 1\) でペナルティが \(q_0 q_1\) の不等式制約を作成する:
```
>>> c = Constraint(q[0] + q[1], le=1, penalty=q[0] * q[1])
```
制約条件式が \(q_0 + q_1 \geq 1\) でペナルティが \(q_0 q_1 - q_0 - q_1\) の不等式制約を作成する:
```
>>> c = Constraint(q[0] + q[1], ge=1, penalty=q[0] * q[1] - q[0] - q[1])
```
制約条件式が \(1 \leq q_0 + q_1 + q_2 \leq 2\) でペナルティが \((q_0 + q_1 + q_2 - 1)(q_0 + q_1 + q_2 - 2)\) の不等式制約を作成する:
```
>>> f = q[0] + q[1] + q[2]
>>> c = Constraint(f, bounds=(1, 2), penalty=(f - 1) * (f - 2))
```