二次割当問題¶
二次割当問題は、以下のような問題です。
\(N\) を正の整数とする。\(N\) 個の工場建設候補地に \(N\) 個の工場を建設することを考える。それぞれの工場は、どの建設候補地に建設されても良い。どの 2 つの工場にも、それらを行き来するトラックが走っており、その輸送量はあらかじめ分かっている。輸送量 x 走行距離の和を最小化するにはどうすればよいか。
応用として、集会の座席表を仲が良い人同士が近くなるように決定することなどが考えられます。
定式化¶
\(N\) 個の工場建設候補地を 土地 \(0\), 土地 \(1\), ..., 土地 \(N-1\) と表記し、\(N\) 個の工場を工場 \(0\), 工場 \(1\), ..., 工場 \(N-1\) と表記します。また土地 \(i\) と土地 \(j\) の距離を \(D_{i, j}\), 工場 \(k\) と工場 \(l\) の間の輸送量を \(F_{k, l}\) とします。
変数¶
\(N \times N\) 個のバイナリ変数 \(q\) を用意し、\(q_{i, k}\) は工場 \(k\) を土地 \(i\) に建設するかどうかを表すことにします。
たとえば、\(q\) が以下のような値をとるとき、土地 \(0\) には 工場 \(3\) が建設されます。
工場 0 |
工場 1 |
工場 2 |
工場 3 |
工場 4 |
|
---|---|---|---|---|---|
土地 0 |
0 |
0 |
0 |
1 |
0 |
土地 1 |
0 |
1 |
0 |
0 |
0 |
土地 2 |
0 |
0 |
0 |
0 |
1 |
土地 3 |
1 |
0 |
0 |
0 |
0 |
土地 4 |
0 |
0 |
1 |
0 |
0 |
制約条件¶
バイナリ変数テーブルのそれぞれの行と列には 1 となる値がちょうど 1 個である必要があります。したがって、各行各列に one-hot 制約をかけます。逆に、これらがみたされていれば、どの土地にどの工場を建設するかが 1 通りに定まります。
目的関数¶
目的関数は、工場と工場の間の輸送量 x 距離の総和です。これは \(q\) を用いて数式で表すと
となります。
数式¶
以上の定式化は、\(N\times N\) 個のバイナリ変数 \(q\) を用いて
と書くことができます。
問題の作成¶
Amplify SDK による定式化を行う前に、問題を作成しておきます。簡単のため工場の数 \(N\) を 10 とします。
import numpy as np
N = 10
土地間の距離を表す行列 \(D\) を作成します。土地はユークリッド平面上にランダムに生成します。距離行列は 2 次元の numpy.ndarray
として作成します。名前は distance
とします。
rng = np.random.default_rng()
x = rng.integers(0, 100, size=(N,))
y = rng.integers(0, 100, size=(N,))
distance = (
(x[:, np.newaxis] - x[np.newaxis, :]) ** 2
+ (y[:, np.newaxis] - y[np.newaxis, :]) ** 2
) ** 0.5
print(distance)
[[ 0. 62.426 12.369 65.947 24.515 49.406 11.18 47.127 58.009 16.031]
[ 62.426 0. 55.154 16.125 85.475 111.83 61.741 65.62 70.114 73.539]
[ 12.369 55.154 0. 61.4 36.056 58.138 21.401 54.589 46.755 18.385]
[ 65.947 16.125 61.4 0. 86.683 114.586 62.482 57.219 84.095 79.322]
[ 24.515 85.475 36.056 86.683 0. 28.844 24.207 50.99 76.694 25.495]
[ 49.406 111.83 58.138 114.586 28.844 0. 52.402 78.435 87.144 41.11 ]
[ 11.18 61.741 21.401 62.482 24.207 52.402 0. 36.139 68.118 26.306]
[ 47.127 65.62 54.589 57.219 50.99 78.435 36.139 0. 99.81 62.37 ]
[ 58.009 70.114 46.755 84.095 76.694 87.144 68.118 99.81 0. 51.264]
[ 16.031 73.539 18.385 79.322 25.495 41.11 26.306 62.37 51.264 0. ]]
工場間の輸送量を表す行列 \(F\) を作成します。2 次元の対称行列をランダムに作成し、名前は flow
とします。
flow = np.zeros((N, N), dtype=int)
for i in range(N):
for j in range(i + 1, N):
flow[i, j] = flow[j, i] = rng.integers(0, 100)
print(flow)
[[ 0 59 2 47 98 90 40 55 40 92]
[59 0 74 20 43 94 7 59 1 24]
[ 2 74 0 91 28 88 84 58 68 76]
[47 20 91 0 94 61 83 14 44 52]
[98 43 28 94 0 45 55 38 32 29]
[90 94 88 61 45 0 14 32 69 4]
[40 7 84 83 55 14 0 97 79 94]
[55 59 58 14 38 32 97 0 95 73]
[40 1 68 44 32 69 79 95 0 26]
[92 24 76 52 29 4 94 73 26 0]]
Amplify SDK による定式化¶
Amplify SDK による定式化を行います。定式化において、どの 2 つのバイナリ変数からなる 2 次項も目的関数に現れるので、Matrix
クラスを使うと効率的な定式化を行うことができます。
変数の作成¶
Matrix
クラスを使用した定式化を行うためには、VariableGenerator
の matrix()
メソッドを使用して変数を発行します。
from amplify import VariableGenerator
gen = VariableGenerator()
matrix = gen.matrix("Binary", N, N) # coefficient matrix
q = matrix.variable_array # variables
q
目的関数の作成¶
上で作成した matrix
は Matrix
クラスのインスタンスであり、
の 3 種類のプロパティを持ちます。
quadratic
は 2 次の項の係数を表す numpy.ndarray
で、その shape
は今回は (N, N, N, N)
となっています。quadratic[i, k, j, l]
は q[i, k] * q[j, l]
の係数に対応します。つまり quadratic
には quadratic[i, k, j, l] = distance[i, j] * flow[k, l]
となるように 4 次元 NumPy 配列を設定する必要があります。
linear
と constant
はそれぞれ線形項の係数と定数項を表しますが、今回使用する目的関数には 2 次の項しか含まれないため、設定しません。
np.einsum("ij,kl->ikjl", distance, flow, out=matrix.quadratic)
制約条件の作成¶
変数の作成
で作成した変数配列 q
の各行各列に one-hot 制約をかけます。
from amplify import one_hot
constraints = one_hot(q, axis=1) + one_hot(q, axis=0)
組合せ最適化モデルの作成¶
目的関数と制約条件を組み合わせて、モデルを作成します。
penalty_weight = np.max(distance) * np.max(flow) * (N - 1)
model = matrix + penalty_weight * constraints
制約条件に penalty_weight
をかけているのは、制約条件に重みをつけるためです。今回使用するソルバーである Amplify AE においては、制約条件に適切な重みを指定しないと、ソルバーが制約条件をみたそうとするよりも目的関数を小さくする方向に動いてしまい、良い解を発見することができなくなってしまいます。詳しくは 制約条件とペナルティ関数 を参照してください。
ソルバークライアントの作成¶
Amplify AE を用いて組合せ最適化を実行するために、ソルバークライアントを作成します。Amplify AE に対応するソルバークライアントクラスは FixstarsClient
クラスです。
from amplify import FixstarsClient
client = FixstarsClient()
Amplify AE の実行に必要な API トークンを設定します。
Tip
ユーザ登録を行うと、評価・検証目的に使える API トークンを無料で入手できます。
client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
ソルバーのタイムアウト値を設定します。
import datetime
client.parameters.timeout = datetime.timedelta(seconds=1)
求解の実行¶
作成した組合せ最適化モデルとソルバークライアントを使用してソルバーの実行を行い、二次計画問題の解を求めます。
from amplify import solve
result = solve(model, client)
目的関数の値は以下のように表示できます。
result.best.objective
239093.17392759357
最適解における変数の値は、NumPy の多次元配列の形式で以下のようにして取得できます。
q_values = q.evaluate(result.best.values)
print(q_values)
[[0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
[0. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]
[0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]
[1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]
[0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]]
結果の確認¶
matplotlib を用いて結果を可視化します。
import itertools
import matplotlib.pyplot as plt
plt.scatter(x, y)
factory_indices = (q_values @ np.arange(N)).astype(int)
for i, j in itertools.combinations(range(N), 2):
plt.plot(
[x[i], x[j]],
[y[i], y[j]],
c="b",
alpha=flow[factory_indices[i], factory_indices[j]] / 100,
)