二次割当問題

二次割当問題は、以下のような問題です。

二次割当問題

\(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\) を用いて数式で表すと

\[ \sum_{q_{i, k} = 1, q_{j, l} = 1} D_{i, j} \ F_{k, l} = \sum_{i, j, k, l} q_{i, k} \ q_{j, l} \ D_{i, j} \ F_{k, l} \]

となります。

数式

以上の定式化は、\(N\times N\) 個のバイナリ変数 \(q\) を用いて

\[\begin{split} \begin{aligned} \text{minimize} \quad &\sum_{i, j, k, l} q_{i, k} \ q_{j, l} \ D_{i, j} \ F_{k, l} \\ \text{subject to} \quad &\sum_k q_{i, k} = 1 \quad \text{for} \quad i \in \{0, 1, \ldots, N - 1\}, \\ &\sum_i q_{i, k} = 1 \quad \text{for} \quad k \in \{0, 1, \ldots, N - 1\}, \\ &q_{i, k} \in \{0, 1\} \quad \text{for} \quad i, k \in \{0, 1, \ldots, N - 1\} \end{aligned} \end{split}\]

と書くことができます。

問題の作成

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.    39.408 22.825 56.851 21.    64.885 39.812 51.264 33.541 42.297]
 [39.408  0.    43.105 71.169 32.062 34.015 28.844 80.156 62.514 76.158]
 [22.825 43.105  0.    34.655 11.045 55.973 26.87  37.054 19.647 34.438]
 [56.851 71.169 34.655  0.    42.72  68.25  44.283 28.636 29.614 39.459]
 [21.    32.062 11.045 42.72   0.    47.508 19.698 48.094 30.594 44.944]
 [64.885 34.015 55.973 68.25  47.508  0.    29.206 88.102 74.061 89.627]
 [39.812 28.844 26.87  44.283 19.698 29.206  0.    59.808 44.944 60.465]
 [51.264 80.156 37.054 28.636 48.094 88.102 59.808  0.    18.248 14.318]
 [33.541 62.514 19.647 29.614 30.594 74.061 44.944 18.248  0.    15.62 ]
 [42.297 76.158 34.438 39.459 44.944 89.627 60.465 14.318 15.62   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 12 31 76 63 73 26 87 29 43]
 [12  0 52 54 34  1 73 17 43 88]
 [31 52  0 21  8  6 83 72  2 65]
 [76 54 21  0 69 10 77 24 74 68]
 [63 34  8 69  0 64 44 43 99 45]
 [73  1  6 10 64  0 96 32 71 14]
 [26 73 83 77 44 96  0 26 77 47]
 [87 17 72 24 43 32 26  0 82 80]
 [29 43  2 74 99 71 77 82  0 30]
 [43 88 65 68 45 14 47 80 30  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

\[\begin{split}\displaystyle \begin{aligned}&\left[\begin{matrix}q_{0,0}& q_{0,1}& q_{0,2}& q_{0,3}& q_{0,4}& q_{0,5}& q_{0,6}& q_{0,7}, q_{0,8}& q_{0,9}\\q_{1,0}& q_{1,1}& q_{1,2}& q_{1,3}& q_{1,4}& q_{1,5}& q_{1,6}& q_{1,7}, q_{1,8}& q_{1,9}\\q_{2,0}& q_{2,1}& q_{2,2}& q_{2,3}& q_{2,4}& q_{2,5}& q_{2,6}& q_{2,7}, q_{2,8}& q_{2,9}\\q_{3,0}& q_{3,1}& q_{3,2}& q_{3,3}& q_{3,4}& q_{3,5}& q_{3,6}& q_{3,7}, q_{3,8}& q_{3,9}\\q_{4,0}& q_{4,1}& q_{4,2}& q_{4,3}& q_{4,4}& q_{4,5}& q_{4,6}& q_{4,7}, q_{4,8}& q_{4,9}\\q_{5,0}& q_{5,1}& q_{5,2}& q_{5,3}& q_{5,4}& q_{5,5}& q_{5,6}& q_{5,7}, q_{5,8}& q_{5,9}\\q_{6,0}& q_{6,1}& q_{6,2}& q_{6,3}& q_{6,4}& q_{6,5}& q_{6,6}& q_{6,7}, q_{6,8}& q_{6,9}\\q_{7,0}& q_{7,1}& q_{7,2}& q_{7,3}& q_{7,4}& q_{7,5}& q_{7,6}& q_{7,7}, q_{7,8}& q_{7,9}\\q_{8,0}& q_{8,1}& q_{8,2}& q_{8,3}& q_{8,4}& q_{8,5}& q_{8,6}& q_{8,7}, q_{8,8}& q_{8,9}\\q_{9,0}& q_{9,1}& q_{9,2}& q_{9,3}& q_{9,4}& q_{9,5}& q_{9,6}& q_{9,7}, q_{9,8}& q_{9,9}\end{matrix}\right]\end{aligned}\end{split}\]

目的関数の作成

上で作成した matrix は Matrix クラスのインスタンスであり、

• quadratic

• linear

• constant

の 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)

組合せ最適化モデルの作成

目的関数と制約条件を組み合わせて、モデルを作成します。

model = matrix + constraints

目的関数と制約条件を足し合わせるだけで、制約条件に重みをつける必要はありません。今回使用する Amplify AE は、制約条件の重みをソルバーの内部で自動的に調整するためです。重みを設定しても Amplify AE には送られず、無視されます。ペナルティ関数を用いて制約条件を扱うほかのソルバーでは、重みの設定が必要です。詳しくは 制約条件とペナルティ関数 を参照してください。

ソルバークライアントの作成

Amplify AE を用いて組合せ最適化を実行するために、ソルバークライアントを作成します。Amplify AE に対応するソルバークライアントクラスは AmplifyAEClient クラスです。

from amplify import AmplifyAEClient

client = AmplifyAEClient()

Amplify AE の実行に必要な API トークンを設定します。

Tip

ユーザ登録を行うと、評価・検証目的に使える API トークンを無料で入手できます。

client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"

ソルバーのタイムアウト値を設定します。

import datetime

client.parameters.time_limit_ms = datetime.timedelta(seconds=1)

求解の実行

作成した組合せ最適化モデルとソルバークライアントを使用してソルバーの実行を行い、二次計画問題の解を求めます。

from amplify import solve

result = solve(model, client)

目的関数の値は以下のように表示できます。

result.best.objective

162442.52084422507

最適解における変数の値は、NumPy の多次元配列の形式で以下のようにして取得できます。

q_values = q.evaluate(result.best.values)

print(q_values)

[[0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]
 [0. 1. 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.]
 [0. 0. 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. 1.]
 [1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 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,
    )