Quadratic Assignment Problem¶
Quadratic assignment problem (QAP) is the following 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.
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.
Formulation¶
The above formulation, with \(N\times N\) binary variables \(q\), can be written as follows.
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. 29.41088234 69.23149572 43.0464865 30.01666204
51.41984053 71.16881339 24.04163056 73.54590403 37.36308338]
[ 29.41088234 0. 78.44743463 71.84705979 24.08318916
61.13100686 100.57832769 5.38516481 82.03657721 58.89821729]
[ 69.23149572 78.44743463 0. 86.55634003 54.5710546
18.02775638 93.64827815 75.66372975 4.47213595 39.40812099]
[ 43.0464865 71.84705979 86.55634003 0. 70.21395873
73. 31.40063694 66.61080993 90.90654542 47.38143096]
[ 30.01666204 24.08318916 54.5710546 70.21395873 0.
37.69615365 94.19129471 22.20360331 58.05170109 41.67733197]
[ 51.41984053 61.13100686 18.02775638 73. 37.69615365
0. 85.0940656 58.05170109 22.20360331 26.07680962]
[ 71.16881339 100.57832769 93.64827815 31.40063694 94.19129471
85.0940656 0. 95.21029356 97.41663102 59.9082632 ]
[ 24.04163056 5.38516481 75.66372975 66.61080993 22.20360331
58.05170109 95.21029356 0. 79.40403012 54.20332093]
[ 73.54590403 82.03657721 4.47213595 90.90654542 58.05170109
22.20360331 97.41663102 79.40403012 0. 43.829214 ]
[ 37.36308338 58.89821729 39.40812099 47.38143096 41.67733197
26.07680962 59.9082632 54.20332093 43.829214 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 13 96 12 75 41 58 23 39 30]
[13 0 72 89 83 87 5 54 61 35]
[96 72 0 38 26 52 19 94 43 44]
[12 89 38 0 4 90 70 56 53 70]
[75 83 26 4 0 99 46 72 47 98]
[41 87 52 90 99 0 74 34 54 7]
[58 5 19 70 46 74 0 33 83 54]
[23 54 94 56 72 34 33 0 47 91]
[39 61 43 53 47 54 83 47 0 45]
[30 35 44 70 98 7 54 91 45 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
Creating the objective function¶
The matrix created above is an instance of the class Matrix, which has the following three properties.
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.
penalty_weight = np.max(distance) * np.max(flow) * (N - 1)
model = matrix + penalty_weight * constraints
The penalty_weight is applied to the constraints to give weight to the constraints. In Amplify AE, the solver used in this example, if you do not specify appropriate weights for the constraints, the solver will search in the direction of making the objective function smaller rather than trying to satisfy the constraints, and you will not be able to find a feasible solution. 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
238294.43673857383
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. 1. 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. 1. 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.]
[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.]]
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,
)