%reload_ext autoreload
%autoreload 2
%matplotlib inline
%load_ext lab_black
割当問題
割当問題(assignment problem)とは，集合 $V=\{1, \ldots, n\}$ および $n \times n$ 行列 $C=[c_{ij}]$ が与えられたとき，
$$ \sum_{i \in V} c_{i, \pi(i)} $$
を最小にする順列 $\pi: V \rightarrow \{1 \ldots, n\}$ を求める問題である．
わかりやすさのために，仕事を資源に割り当てる問題として解釈しておこう． $n$ 個の仕事があり，それを $n$ 個の資源（作業員と考えても良い）に割り当てることを考える． 作業と資源には相性があり，仕事 $i$ を資源 $j$ に割り当てた際の費用を $c_{ij}$ とする． 1つの作業員は高々1つの仕事しか処理できず，1つの仕事には必ず1つの資源を割り当てるとき，総費用を最小化する問題が割当問題である．
仕事 $i$ を資源 $j$ に割り当てるとき $1$，それ以外のとき $0$ となる $0$-$1$ 変数を使うと，割当問題は以下のように整数最適化問題として定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{i,j \in V} c_{ij} x_{ij} & \\ s.t. & \sum_{j \in V} x_{ij} = 1 & \forall i \in V \\ & \sum_{i \in V} x_{ij} = 1 & \forall j \in V \\ & x_{ij} \in \{0,1\} & \forall i,j \in V \end{array} $$この問題は，変数を実数に緩和しても，最適解においては整数になることが知られている．これを完全単模性(total modularity)とよぶ． したがって， $x_{ij} \geq 0$ に緩和して，線形最適化ソルバーで求解すれば，最適解を得ることができる． また，問題を変形することによって，最小費用流問題に帰着できるので，networkXを用いても解くことができる．最小費用流問題の解法には，ネットワーク単体法と容量スケーリング法があるが，前者の方が高速である． また，networkXの実装はPythonであるので，Gurobiなどの商用の数理最適化ソルバーを用いた方が高速である．なお，この問題は極度に退化しているので，退化対策が十分になされていない単体法の実装では，時間がかかることがある．
SciPyにHungary法の実装があるので，これを使う方が高速である．
networkXの最小費用流を用いても解けるが，network_simplexは費用が浮動小数点数の場合には非常に遅くなるので，数値を整数に丸めてから実行する必要がある． capacity_scalingは，費用が浮動小数点数でも性能に差はないが，network_simplexより遅い．
n = 1000
cost = np.random.randint(100, 1000, size=(n, n))
%%time
row_ind, col_ind = linear_sum_assignment(cost)
print(cost[row_ind, col_ind].sum())
G = nx.DiGraph()
n = len(cost)
for i in range(n):
G.add_node(i, demand=-1)
G.add_node(n + i, demand=1)
G.add_weighted_edges_from([(i, n + j, cost[i, j]) for i in range(n) for j in range(n)])
%%time
val, flow = nx.algorithms.flow.network_simplex(G)
print(val)
%%time
val, flowDict = nx.capacity_scaling(G)
print(val)
V = list(range(n))
model = Model("ap")
x = {}
for i in V:
for j in V:
x[i, j] = model.addVar(vtype="C", 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(cost[i, j] * x[i, j] for i in V for j in V), GRB.MINIMIZE)
%%time
model.optimize()
print(model.ObjVal)
ボトルネック割当問題
ボトルネック割当問題 (bottleneck assignment problem) とは，集合 $V=\{1, \cdots, n\}$ および $n \times n$ 行列 $C=[c_{ij}]$ が与えられたとき，
$$ \max_{i \in V} c_{i \pi(i)} $$を最小にする順列 $\pi: V \rightarrow \{1 \cdots, n\}$ を求める問題である．
仕事を資源に割り当てる例においては，行列 $C$ を処理時間としたとき，最も遅く終了する時間を最小化することに相当する． 仕事 $i$ が資源 $j$ に割り当てられるとき $1$，それ以外のとき $0$ となる $0$-$1$ 変数 $x_{ij}$ を用いると，ボトルネック割当問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & z & \\ s.t. & \sum_{i \in V} c_{ij} x_{ij} \leq z & \forall j \in V \\ & \sum_{i \in V} x_{ij} = 1 & \forall j \in V \\ & \sum_{j \in V} x_{ij} = 1 & \forall i \in V \\ & x_{ij} \in \{0,1\} & \forall i,j \in V \end{array} $$V = list(range(n))
model = Model("bap")
x = {}
for i in V:
for j in V:
x[i, j] = model.addVar(vtype="B", name=f"x[{i},{j}]")
z = model.addVar(vtype="C", name="z")
model.update()
for i in V:
model.addConstr(quicksum(cost[i, j] * x[i, j] for j in V) <= z)
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(z, GRB.MINIMIZE)
%%time
model.optimize()
print("obj=", z.X)
%%time
LB = max(cost.min(axis=1).max(), cost.min(axis=0).max() )
row_ind, col_ind = linear_sum_assignment(cost)
UB =cost[row_ind, col_ind].max()
for t in range(LB, UB):
c = cost.flatten()
c = np.where(c > t, np.inf, c)
c.shape=(n,n)
try:
row_ind, col_ind = linear_sum_assignment(c)
break
except:
continue
print(cost[row_ind, col_ind].max())
一般化割当問題
$n$ 個の仕事があり，それを $m$ 個の資源（作業員と考えても良い）に割り当てることを考える． 作業と資源には相性があり，仕事 $j$ を資源 $i$ に割り当てた際の費用を $c_{ij}$ とする． 資源 $i$ には使用可能量の上限（容量） $b_i$ が定義されており，仕事 $j$ を資源 $i$ に割り当てるときには，資源 $a_{ij}$ を使用とするものする． 1つの仕事には必ず1つの資源を割り当てるとき，資源に割り当てられた仕事の資源使用量の合計が資源の容量を超えないという条件のもとで，総費用を最小化する問題が一般化割当問題(generalized assignment problem)である．
仕事 $j$ を資源 $i$ に割り当てるとき $1$，それ以外のとき $0$ となる $0$-$1$ 変数 $x_{ij}$ を用いると，一般化割当問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{i,j} c_{ij} x_{ij} & \\ s.t. & \sum_{j} a_{ij} x_{ij} \leq b_i & \forall i \\ & \sum_{i} x_{ij} = 1 & \forall j \\ & x_{ij} \in \{0,1\} & \forall i,j \end{array} $$一般化割当問題は，ナップサック問題を含むので，$NP$-困難である．
以下のサイトからダウンロードでいるベンチマーク問題例を用いる．
https://www-or.amp.i.kyoto-u.ac.jp/members/yagiura/gap/
以下に，データ読み込みのコードを示す．
folder = "../data/gap/"
f = open(folder + "a05100")
data = f.readlines()
f.close()
m, n = list(map(int, data[0].split()))
cost, a, b = {}, {}, {}
data_ = []
for row in data[1:]:
data_.extend(list(map(int, row.split())))
count = 0
for i in range(m):
for j in range(n):
cost[i, j] = data_[count]
count += 1
for i in range(m):
for j in range(n):
a[i, j] = data_[count]
count += 1
for i in range(m):
b[i] = data_[count]
count += 1
メタヒューリスティクス
以下の論文に掲載されているメタヒューリスティクスをPythonから呼び出して使用する．
M. Yagiura, T. Ibaraki and F. Glover, An Ejection Chain Approach for
the Generalized Assignment Problem, INFORMS Journal on Computing, 16
(2004) 133-151. M. Yagiura, T. Ibaraki and F. Glover, A Path Relinking Approach with
Ejection Chains for the Generalized Assignment Problem, European
Journal of Operational Research, 169 (2006) 548-569.このコード自体はオープンソースではないが，研究用なら筆者に連絡すれば利用可能である． 詳細は，付録1の OptGAP を参照されたい．
file_name = "a05100"
timelimit = 1
sol = gapmeta(folder + file_name, timelimit)
print("solution=", sol)
I = list(range(m))
J = list(range(n))
model = Model("gap")
x = {}
for i in I:
for j in J:
x[i, j] = model.addVar(vtype="B", name=f"x[{i},{j}]")
model.update()
for i in I:
model.addConstr(quicksum(a[i, j] * x[i, j] for j in J) <= b[i])
for j in J:
model.addConstr(quicksum(x[i, j] for i in I) == 1)
model.setObjective(quicksum(cost[i, j] * x[i, j] for i in I for j in J), GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)