二次割当問題

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

二次割当問題

\(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{align} \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{align} \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.    84.054 76.118 43.139 35.903 29.682 90.25  96.799 72.443 23.087]
 [84.054  0.    18.358 63.812 49.254 61.984 36.    30.806 13.    62.801]
 [76.118 18.358  0.    48.466 45.177 50.09  21.932 48.703 21.213 53.339]
 [43.139 63.812 48.466  0.    39.319 15.811 54.037 88.414 57.315 25.298]
 [35.903 49.254 45.177 39.319  0.    26.    63.891 61.327 37.    19.849]
 [29.682 61.984 50.09  15.811 26.     0.    61.4   82.201 52.773  9.487]
 [90.25  36.    21.932 54.037 63.891 61.4    0.    66.37  42.72  67.231]
 [96.799 30.806 48.703 88.414 61.327 82.201 66.37   0.    31.401 79.756]
 [72.443 13.    21.213 57.315 37.    52.773 42.72  31.401  0.    52.202]
 [23.087 62.801 53.339 25.298 19.849  9.487 67.231 79.756 52.202  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 62 34 52  1  2 97 34 68 76]
 [62  0 82 23 77 11 40 20 40 54]
 [34 82  0  2 97 85 19 71 11 60]
 [52 23  2  0 56 35 72 30 30 94]
 [ 1 77 97 56  0 96  9  1 86 64]
 [ 2 11 85 35 96  0 19  7 17 82]
 [97 40 19 72  9 19  0 96  1 90]
 [34 20 71 30  1  7 96  0 23 80]
 [68 40 11 30 86 17  1 23  0 29]
 [76 54 60 94 64 82 90 80 29  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)

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

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

penalty_weight = np.max(distance) * np.max(flow) * (N - 1)
model = matrix + penalty_weight * constraints

制約条件に penalty_weight をかけているのは、制約条件に重みをつけるためです。今回使用するソルバーである 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

165115.3786323302

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

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

print(q_values)

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