2次割当問題とは
2次割当問題(quadratic assignment problem)とは，集合 $V=\{1, \cdots, n\}$ および2つの $n \times n$ 行列 $F=[f_{ij}]$，$D=[d_{k \ell}]$ が与えられたとき，
$$ \sum_{i \in V} \sum_{j \in V} f_{ij} d_{\pi(i) \pi(j)} $$
を最小にする順列 $\pi: V \rightarrow \{1 \ldots, n\}$ を求める問題である．
この問題は,施設の配置を決定する応用から生まれた．$n$ 個の施設があり，それを $n$ 箇所の地点に配置することを 考える．施設 $i,j$ 間には物の移動量 $f_{ij}$ があり， 地点 $k,\ell$ 間を移動するには距離 $d_{k \ell}$ がかかるものとする． 問題の目的は，物の総移動距離を最小にするように， 各地点に1つずつ施設を配置することである． 順列 $\pi$ は施設 $i$ を地点 $\pi(j)$ に配置することを表す．
2次割当問題にはQAPLIBとよばれるベンチマーク問題集がある．以下のサイトから入手できる．
folder = "../data/qap/"
df = pd.read_csv(folder + "qaplib.csv")
df
定式化
施設 $i$ が地点 $k$ に配置されるとき $1$，それ以外のとき $0$ となる $0$-$1$ 変数 $x_{ik}$ を用いると，2次割当問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{i,j \in V} \sum_{k, \ell \in V} f_{ij} d_{k \ell} x_{ik} x_{j \ell} & \\ s.t. & \sum_{i \in V} x_{ik} = 1 & \forall k \in V \\ & \sum_{k \in V} x_{ik} = 1 & \forall i \in V \\ & x_{ik} \in \{0,1\} & \forall i,k \in V \end{array} $$市販の数理最適化ソルバーは，通常，上のような（下に凸でない）2次の関数を含んだ最小化問題には対応していなかったが，最近ではGurobiがこれに対応した解法を取り入れた． 小規模で疎な行列の問題例なら，比較的短時間で解ける場合もある．
整数線形最適化に帰着させるためには，様々な方法が提案されているが，ここでは2つの定式化を紹介しよう．
はじめの定式化は，施設 $i$ を地点 $k$ に配置したときの，施設 $i$ と他の施設の間の費用の合計を表す実数変数 $w_{ik}$ を追加したものであり， $O(n^2)$ 個の変数を用いたコンパクトなものである．
施設 $i$ を地点 $k$ に配置したときの費用の上界となるパラメータ $M_{ik}$ を導入する． $$ M_{ik} =\sum_{j,\ell \in V} f_{ij} d_{k \ell} $$
線形整数最適化による定式化は以下のように書ける．
$$ \begin{array}{l l l} minimize & \sum_{i,k \in V} w_{ik} & \\ s.t. & \sum_{i \in V} x_{ik} = 1 & \forall k \in V \\ & \sum_{k \in V} x_{ik} = 1 & \forall i \in V \\ & M_{ik}(x_{ik}-1) +\sum_{j,\ell \in V} f_{ij} d_{k \ell} x_{j \ell} \leq w_{ik} & \forall i,k \in V \\ & x_{ik} \in \{0,1\} & \forall i,k \in V \\ & w_{ik} \geq 0 & \forall i,k \in V \end{array} $$次の定式化は，施設 $i$ が地点 $k$ に配置され，かつ施設 $j$ が地点 $\ell$に配置されるとき $1$，それ以外のとき $0$ となる4添え字の $0$-$1$ 変数 $y_{ikj\ell}$ を用いたものである． 変数の数は $O(n^4)$ 個と多いが，得られる下界は強くなる． 線形整数最適化による定式化は以下のように書ける．
$$ \begin{array}{l l l} minimize & \sum_{i,j \in V, i \neq j} \sum_{k, \ell \in V, k \neq \ell} f_{ij} d_{k \ell} y_{ikj \ell} & \\ s.t. & \sum_{i \in V} x_{ik} = 1 & \forall k \in V \\ & \sum_{k \in V} x_{ik} = 1 & \forall i \in V \\ & \sum_{i \in V, i \neq j} y_{ikj\ell} =x_{j\ell} & \forall k,j,\ell \in V, k \neq \ell \\ & \sum_{k \in V, k \neq \ell} y_{ikj\ell} =x_{j\ell} & \forall i,j,\ell \in V, i \neq j \\ & y_{ikj\ell}=y_{j\ell ik} & \forall i,j,k,\ell \in V, i \neq j, k \neq \ell \\ & x_{ik} \in \{0,1\} & \forall i,k \in V \\ & y_{ikj\ell} \in \{0,1\} & \forall i,k,j,\ell \in V, i \neq j, k \neq \ell \end{array} $$以下ではベンチマーク問題例を読み込み，$O(n^2)$ 個の変数を用いた定式化を用いて数理最適化ソルバーで求解する． ただし，この問題は難しいので，計算を途中で打ち切り，得られた最良解で評価する．
def read_qap(filename):
"""Read data for a QAP problem from file in QAPLIB format."""
try:
if len(filename) > 3 and filename[-3:] == ".gz": # file compressed with gzip
import gzip
f = gzip.open(filename, "rb")
else: # usual, uncompressed file
f = open(filename)
except IOError:
print("could not open file", filename)
return None
data = f.read()
f.close()
try:
pass
data = data.split()
n = int(data.pop(0))
f = {} # for n times n flow matrix
d = {} # for n times n distance matrix
for i in range(n):
for j in range(n):
f[i, j] = int(data.pop(0))
for i in range(n):
for j in range(n):
d[i, j] = int(data.pop(0))
except IndexError:
print("inconsistent data on QAP file", filename)
exit(-1)
return n, f, d
n, f, d = read_qap(folder +"tai20a.dat")
V = list(range(n))
model = Model("qap")
x = {}
for i in V:
for j in V:
x[i, j] = model.addVar(vtype="B", name=f"x[{i},{j}]")
model.update()
for j in V:
model.addConstr(quicksum(x[i, j] for i in V) == 1)
for i in V:
model.addConstr(quicksum(x[i, j] for j in V) == 1)
model.setObjective(
quicksum(
f[i, j] * d[k, ell] * x[i, k] * x[j, ell]
for i in V
for j in V
for k in V
for ell in V
),
GRB.MINIMIZE,
)
model.Params.NoRelHeurTime = 600.
#model.Params.PreQLinearize = 1
model.optimize()
n, f, d = read_qap(folder +"tai20a.dat")
V = list(range(n))
M ={}
for i in V:
for k in V:
sum_ = 0
for j in V:
for ell in V:
sum_ += f[i,j]*d[k,ell]
M[i,k] = sum_
model = Model("qap-linear")
x, w = {}, {}
for i in V:
for j in V:
w[i, j] = model.addVar(vtype="C", name=f"w[{i},{j}]")
x[i, j] = model.addVar(vtype="B", name=f"x[{i},{j}]")
model.update()
for j in V:
model.addConstr(quicksum(x[i, j] for i in V) == 1)
for i in V:
model.addConstr(quicksum(x[i, j] for j in V) == 1)
for i in V:
for k in V:
model.addConstr(M[i,k]*(x[i,k]-1) +
quicksum(f[i,j]*d[k,ell]*x[j,ell] for j in V for ell in V) <= w[i,k])
model.setObjective(
quicksum(
w[i, k]
for i in V
for k in V
),
GRB.MINIMIZE,
)
model.Params.NoRelHeurTime = 600.
#model.Params.PreQLinearize = 1
model.optimize()
タブーサーチ
2次割当問題に対しては，解を表す順列の2つの要素を交換することによって 自然な近傍が定義できる．これは，2つの施設の配置場所を交換することを表す． これを交換近傍とよぶ．
交換近傍と補助行列による目的関数値の差分の高速計算を用いて，タブーサーチを設計する． タブーサーチ対しては，何を禁断リストに入れるのかを決める必要がある． 定式化では，施設 $i$ を地点 $k$ に配置するときに $1$ になる変数 $x_{ik}$ を用いていたことから， 施設 $i$ を地点 $k$ から他の地点に移動したとき，再び地点 $k$ に戻ることを一定の反復の間禁止することにする．
集中化と多様化を入れたタブーサーチで，数理最適化で解けなかった問題例を解いてみる．
n, f, d = read_qap(folder +"tai60a.dat")
print("starting tabu search")
tabulen = n
max_iterations = 30000
pi, cost = tabu_search(n, f, d, max_iterations, tabulen, report=print)
print("final solution: z =", cost)
print(pi)
SciPyの近似解法
SciPyにも2次割当問題の近似解法が含まれている． 引数 method で解法を選ぶことができる． 引数には "faq" (Fast Approximate QAP Algorithm)と"2opt"が使える． どちらも高速ではあるが，解はあまり良くない．
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.quadratic_assignment.html
n, f, d = read_qap(folder + "wil50.dat")
n, f, d = read_qap(folder +"tai20a.dat")
F = np.zeros((n, n))
for (i, j) in f:
F[i, j] = f[i, j]
D = np.zeros((n, n))
for (i, j) in d:
D[i, j] = d[i, j]
quadratic_assignment(A=F, B=D, method="faq")
quadratic_assignment(A=F, B=D, method="2opt")
線形順序付け問題
2次割当問題の特殊形に線形順序付け問題 (linear ordering problem) がある．
線形順序付け問題とは，集合 $V=\{1, \ldots, n\}$ および $n \times n$ 行列 $C=[c_{ij}]$ が与えられたとき，
$$ \sum_{i < j } c_{\pi(i) \pi(j)} $$を最小にする順列 $\pi: V \rightarrow \{1 \ldots, n\}$ を求める問題である．
定式化
$\pi (i) < \pi (j)$ のとき $1$ になる $0$-$1$ 変数 $x_{ij}$ を用いると，線形順序付け問題は以下のように定式化できる．
$$ \begin{array}{l l l} maximize & \sum_{i,j \in V, i \neq j } c_{ij} x_{ij} & \\ s.t. & x_{ij} + x_{ji} = 1 & \forall i,j \in V, i < j \\ & x_{ij} + x_{jk} + x_{ki} \leq 2 & \forall i < j, i < k, j \neq k \\ & x_{ij} \in \{0,1\} & \forall i,j \in V, i \neq j \end{array} $$線形順序付け問題のベンチマーク問題例は，以下から入手できる．
folder = "../data/lop/"
fn = "N-be75np"
f = open(folder + fn)
data = f.readlines()
f.close()
n = int(data[0])
print(n)
c = []
for row in data[1:]:
c.append(list(map(int, row.split())))
cost = np.array(c)
cost
sol = "7 48 49 33 38 50 11 28 25 27 41 45 37 23 22 20 24 26 42 43 1 5 21 39 34 19 15 18 16 17 14 32 8 9 10 35 36 30 46 29 31 13 12 3 44 2 40 4 6 47"
sol_list = list(map(int, sol.split()))
total = 0
for i in range(n - 1):
for j in range(i + 1, n):
total += cost[sol_list[i] - 1, sol_list[j] - 1]
total
V = list(range(n))
model = Model("lop")
x = {}
for i in V:
for j in V:
if i != j:
x[i, j] = model.addVar(vtype="B", name=f"x[{i},{j}]")
model.update()
for i in range(n - 1):
for j in range(i + 1, n):
model.addConstr(x[i, j] + x[j, i] == 1)
for i in range(n - 1):
for j in range(i + 1, n):
for k in range(i + 1, n):
if j != k:
model.addConstr(x[i, j] + x[j, k] + x[k, i] <= 2)
model.setObjective(
quicksum(cost[i, j] * x[i, j] for i in V for j in V if i != j), GRB.MAXIMIZE
)
model.optimize()
print(model.ObjVal)
D = nx.DiGraph()
for (i, j) in x:
if x[i, j].X > 0:
D.add_edge(i, j)
perm = []
for i in nx.topological_sort(D):
perm.append(i)
print(perm)