二次割当問題

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

二次割当問題

\(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.     74.431  17.117  34.132  20.248  37.216  36.249  91.81   64.885  15.232]
 [ 74.431   0.     57.315  87.71   57.009  83.217  41.593  42.107  28.601  82.97 ]
 [ 17.117  57.315   0.     41.183   9.22   41.11   20.616  76.837  49.406  28.018]
 [ 34.132  87.71   41.183   0.     34.828   7.616  60.803 116.     88.459  20.224]
 [ 20.248  57.009   9.22   34.828   0.     33.541  26.077  81.32   53.74   25.962]
 [ 37.216  83.217  41.11    7.616  33.541   0.     59.54  113.569  86.215  25.08 ]
 [ 36.249  41.593  20.616  60.803  26.077  59.54    0.     56.223  28.844  48.6  ]
 [ 91.81   42.107  76.837 116.     81.32  113.569  56.223   0.     27.586 104.809]
 [ 64.885  28.601  49.406  88.459  53.74   86.215  28.844  27.586   0.     77.421]
 [ 15.232  82.97   28.018  20.224  25.962  25.08   48.6   104.809  77.421   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 56 51 81 14 32  5 19 42 80]
 [56  0 75 15 93 73 48  3 72 47]
 [51 75  0 16 47 18 21 17 94 55]
 [81 15 16  0 87  7 76 11 79 95]
 [14 93 47 87  0 24 66 62 95 22]
 [32 73 18  7 24  0 63 67 91 48]
 [ 5 48 21 76 66 63  0 67 96 49]
 [19  3 17 11 62 67 67  0 10 38]
 [42 72 94 79 95 91 96 10  0 57]
 [80 47 55 95 22 48 49 38 57  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

190260.64186051858

最適解における変数の値は、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. 0. 0. 0. 0. 0. 0. 1. 0.]
 [1. 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. 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. 0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]]

結果の確認

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