Quadratic Assignment Problem

Quadratic assignment problem (QAP) is the following problem.

Quadratic assignment problem

Let \(N\) be a positive integer. Consider \(N\) factories to be built on \(N\) candidate sites. Each factory can be built on any of the candidate sites. Every two factories have trucks traveling to and from them, and their transportation volumes are known in advance. How can we minimize the sum of the amount transported x the distance traveled?

An application could be to determine the seating chart for a meeting so that people close to each other have seats closely.

Formulation

Let \(N\) potential factory locations be denoted by land \(0\), land \(1\), … , and \(N\) factories are denoted as factories \(0\), factories \(1\), …, factories \(N-1\). Also let \(D_{i, j}\) denote the distance between land \(i\) and land \(j\), and \(F_{k, l}\) denote the transport volume between factory \(k\) and factory \(l\).

Variables

With \(N \times N\) binary variables \(q\), let \(q_{i, k}\) represent whether factory \(k\) is to be built on land \(i\).

For example, factory \(3\) will be built on land \(0\) if \(q\) has the following value.

Binary variable table

	factory 0	factory 1	factory 2	factory 3	factory 4
land 0	0	0	0	1	0
land 1	0	1	0	0	0
land 2	0	0	0	0	1
land 3	1	0	0	0	0
land 4	0	0	1	0	0

Constraints

Each row and column of the binary variable table must have exactly one variable that is 1, so we place a one-hot constraint on each row and column. Conversely, if these are satisfied, then there is only one way to determine which factory to build on which land.

Objective function

The objective function is the sum of transport volume x distance between factories. This can be expressed in the equation using \(q\) as follows.

\[ \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} \]

Formulation

The above formulation, with \(N\times N\) binary variables \(q\), can be written as follows.

\[\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}\]

Problem setting

Before formulating with the Amplify SDK, we will create a problem. For simplicity, let the number of factories \(N=10\).

import numpy as np

N = 10

Next, we create a distance matrix \(D\) representing the distances between lands. The lands are randomly generated on the Euclidean plane. The distance matrix distance is created as a two-dimensional numpy.ndarray.

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.      5.385 115.278  37.643  55.714  63.135  59.405  19.209  65.123  36.892]
 [  5.385   0.    113.172  38.471  56.045  62.65   58.31   14.142  67.602  36.249]
 [115.278 113.172   0.     84.38   68.942  57.697  57.689 114.403  88.284  79.812]
 [ 37.643  38.471  84.38    0.     18.358  27.586  26.833  48.166  32.65   10.296]
 [ 55.714  56.045  68.942  18.358   0.     11.402  15.524  64.195  27.295  21.   ]
 [ 63.135  62.65   57.697  27.586  11.402   0.      8.062  68.884  36.401  26.401]
 [ 59.405  58.31   57.689  26.833  15.524   8.062   0.     63.246  42.544  22.672]
 [ 19.209  14.142 114.403  48.166  64.195  68.884  63.246   0.     79.31   43.29 ]
 [ 65.123  67.602  88.284  32.65   27.295  36.401  42.544  79.31    0.     41.617]
 [ 36.892  36.249  79.812  10.296  21.     26.401  22.672  43.29   41.617   0.   ]]

Also, we create a matrix \(F\) representing the amount of transport between factories, a random symmetric matrix of dimension 2, named 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 24 24 54 33 78  3  0 98 49]
 [24  0 55 21 13 20 50 98 68 94]
 [24 55  0 32 19 55 79 11  2 83]
 [54 21 32  0 64 61 79  6 30 68]
 [33 13 19 64  0  2 86 13  9 43]
 [78 20 55 61  2  0  8 69 63 60]
 [ 3 50 79 79 86  8  0 41 39 55]
 [ 0 98 11  6 13 69 41  0  2 43]
 [98 68  2 30  9 63 39  2  0 96]
 [49 94 83 68 43 60 55 43 96  0]]

Formulation with the Amplify SDK

In the formulation, we can use the Matrix class for efficient formulation, since a quadratic term consisting of any two binary variables can appear in the objective function.

Creating variables

To formulate using the Matrix class, VariableGenerator’s matrix() method to issue variables.

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}\]

Creating the objective function

The matrix created above is an instance of the class Matrix, which has the following three properties.

• quadratic

• linear

• constant

quadratic is numpy.ndarray representing the coefficients of the second order terms, and its shape is (N, N, N, N) this time. quadratic[i, k, j, l] corresponds to the coefficients of q[i, k] * q[j, l]. That is, quadratic must be set to a 4-dimensional NumPy array such that quadratic[i, k, j, l] = distance[i, j] * flow[k, l]

linear and constant represent the coefficient and constant terms of the linear term, respectively, but since the objective function used in this problem contains only second order terms, we will not set them.

np.einsum("ij,kl->ikjl", distance, flow, out=matrix.quadratic)

Creating constraints

Impose a one-hot constraint on each row and column of the variable array q created in Creating variables.

from amplify import one_hot

constraints = one_hot(q, axis=1) + one_hot(q, axis=0)

Creating a combinatorial optimization model

Let’s combine the objective function and constraints to create a model.

model = matrix + constraints

You add the objective function and the constraints together, and you give no weight to the constraints. Amplify AE, the solver of this example, adjusts the weight of each constraint by itself. It does not receive a weight that you set, and it ignores that weight. A solver that uses a penalty function for a constraint needs the weight. See Constraints and Penalty Functions for details.

Creating a solver client

Now, we will create a solver client to perform combinatorial optimization using Amplify AE. The solver client class corresponding to Amplify AE is AmplifyAEClient class.

from amplify import AmplifyAEClient

client = AmplifyAEClient()

We also need to set the API token required to run Amplify AE.

Tip

After user registration, you can obtain a free API token that can be used for evaluation and validation purposes.

client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"

We will set the solver’s timeout.

import datetime

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

Executing the solver

Finally, we will execute the solver using the created combinatorial optimization model and the solver client to find the solution to the quadratic programming problem.

from amplify import solve

result = solve(model, client)

The objective function value based on the best solution is shown below.

result.best.objective

150491.03255839928

The values of the variables in the optimal solution can be obtained in the form of a NumPy multidimensional array as follows.

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

print(q_values)

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

Checking the results

We will visualize the results using 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,
    )