Quadratic Assignment Problem¶
Quadratic assignment problem (QAP) is the following problem.
Let
An application could be to determine the seating chart for a meeting so that people close to each other have seats closely.
Formulation¶
Let
Variables¶
With
For example, factory
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
Formulation¶
The above formulation, with
Problem setting¶
Before formulating with the Amplify SDK, we will create a problem. For simplicity, let the number of factories
import numpy as np
N = 10
Next, we create a 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. 22. 80.056 74.25 13.601 82.377 40.706 18.028 39.812 13. ]
[ 22. 0. 91.831 86.052 13.601 101.518 61.033 4.123 53.075 34.366]
[ 80.056 91.831 0. 5.831 78.816 43.417 53.009 88.6 40.249 80.362]
[ 74.25 86.052 5.831 0. 73.007 43.046 48.166 82.801 34.438 74.726]
[ 13.601 13.601 78.816 73.007 0. 88.142 48.27 9.899 39.598 26.42 ]
[ 82.377 101.518 43.417 43.046 88.142 0. 42.012 97.494 52.773 75.822]
[ 40.706 61.033 53.009 48.166 48.27 42.012 0. 56.921 25.495 33.941]
[ 18.028 4.123 88.6 82.801 9.899 97.494 56.921 0. 49.497 30.594]
[ 39.812 53.075 40.249 34.438 39.598 52.773 25.495 49.497 0. 41.304]
[ 13. 34.366 80.362 74.726 26.42 75.822 33.941 30.594 41.304 0. ]]
Also, we create a matrix 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 76 32 7 97 32 89 68 87 99]
[76 0 47 41 90 78 11 75 64 0]
[32 47 0 57 2 27 2 14 71 90]
[ 7 41 57 0 77 11 87 29 79 82]
[97 90 2 77 0 59 68 19 73 98]
[32 78 27 11 59 0 28 83 52 92]
[89 11 2 87 68 28 0 93 98 41]
[68 75 14 29 19 83 93 0 4 65]
[87 64 71 79 73 52 98 4 0 87]
[99 0 90 82 98 92 41 65 87 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 FixstarsClient
class.
from amplify import FixstarsClient
client = FixstarsClient()
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.timeout = 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
221348.4238258041
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. 0. 0. 0. 1. 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. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]
[0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]
[1. 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. 0.]
[0. 0. 0. 0. 0. 0. 1. 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,
)
