k-センター問題に対する定式化と解法
import random
import math
from gurobipy import Model, quicksum, GRB
# from mypulp import Model, quicksum, GRB
import networkx as nx
import matplotlib.pyplot as plt
標準定式化
$k$-メディアン問題と同様に，顧客 $i$ から施設 $j$ への距離を $c_{ij}$ とし，以下の変数を導入する．
$$ x_{ij}= \left\{ \begin{array}{ll} 1 & 顧客 i の需要が施設 j によって満たされるとき \\ 0 & それ以外のとき \end{array} \right. $$$$ y_j = \left\{ \begin{array}{ll} 1 & 施設 j を開設するとき \\ 0 & それ以外のとき \end{array} \right. $$最も遠い施設でサービスを受ける顧客の移動費用を表す連続変数 $z$ を導入する．
上の記号および変数を用いると，$k$-センター問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & z & \\ s.t. & \sum_{j \in J} x_{ij} =1 & \forall i \in I \\ & \sum_{j \in J} c_{ij} x_{ij} \leq z & \forall i \in I \\ & \sum_{j \in J} y_{j} = k & \\ & x_{ij} \leq y_j & \forall i \in I, j \in J \\ & x_{ij} \in \{ 0,1 \} & \forall i \in I, j \in J \\ & y_j \in \{ 0,1 \} & \forall j \in J \end{array} $$def distance(x1, y1, x2, y2):
return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
def make_data(n):
I = range(n)
J = range(n)
x = [random.random() for i in range(n)]
y = [random.random() for i in range(n)]
c = {}
for i in I:
for j in J:
c[i, j] = distance(x[i], y[i], x[j], y[j])
return I, J, c, x, y
def kcenter(I, J, c, k):
model = Model("k-center")
z = model.addVar(vtype="C")
x, y = {}, {}
for j in J:
y[j] = model.addVar(vtype="B")
for i in I:
x[i,j] = model.addVar(vtype="B")
model.update()
for i in I:
model.addConstr(quicksum(x[i,j] for j in J) == 1)
model.addConstr(quicksum(c[i,j]*x[i,j] for j in J) <= z)
for j in J:
model.addConstr(x[i,j] <= y[j])
model.addConstr(quicksum(y[j] for j in J) == k)
model.setObjective(z, GRB.MINIMIZE)
model.update()
model.__data = x,y
return model
random.seed(67)
n = 30
I, J, c, x_pos, y_pos = make_data(n)
k = 5
model = kcenter(I, J, c, k)
model.optimize()
x, y = model.__data
facilities = [j for j in y if y[j].X > 0.5]
edges = [(i, j) for (i, j) in x if x[i, j].X > 0.5]
def draw_center(facilities, edges, x_pos, y_pos):
G = nx.Graph()
facilities = set(facilities)
unused = set(j for j in J if j not in facilities)
client = set(i for i in I if i not in facilities and i not in unused)
G.add_nodes_from(facilities)
G.add_nodes_from(client)
G.add_nodes_from(unused)
for (i, j) in edges:
G.add_edge(i, j)
position = {}
for i in range(len(x_pos)):
position[i] = (x_pos[i], y_pos[i])
nx.draw(G, position, with_labels=False, node_color="b", nodelist=facilities)
nx.draw(G, position, with_labels=False, node_color="c", nodelist=unused, node_size=50)
nx.draw(G, position, with_labels=False, node_color="g", nodelist=client, node_size=50)
draw_center(facilities, edges, x_pos, y_pos)
2分探索を用いた解法
標準定式化はmin-max型の目的関数をもつので大規模問題例には不向きである． ここでは，2分探索を用いた解法を考える．
顧客から施設までの距離が $\theta$ 以下である枝だけから構成されるグラフ $G_{\theta}=(V,E_{\theta})$ を考える． 点の部分集合 $S(\subseteq V)$ が与えられたとき， すべての点 $i(\in V)$ が $S$ 内の少なくとも1つの点に隣接するとき $S$ を点被覆(vertex cover)と呼ぶ．
グラフ $G_{\theta}$ 上で $|S|=k$ の被覆が存在するなら，$k$-センター問題の最適値が $\theta$ 以下であることが言える．
施設 $j$ を開設するとき $1$ となる $0$-$1$ 変数 $y_j$ は， ここでは点の部分集合 $S$ に含まれるとき $1$ となる変数と見なすことができる． さらに以下の変数を導入する． $$ z_i = \left\{ \begin{array}{ll} 1 & 点 i が S 内の点に隣接しない(被覆されない） \\ 0 & それ以外のとき \end{array} \right. $$
グラフ $G_{\theta}$ の隣接行列（点 $i,j$ が隣接しているとき $1$， それ以外のとき $0$ の要素をもつ行列）を $[a_{ij}]$ としたとき， $|S|=k$ の被覆が存在するか否かを判定する問題は， 以下の被覆立地問題(covering location problem)として定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{i \in I} z_{i} & \\ s.t. & \sum_{j \in J} a_{ij} y_{j} +z_{i} =1 & \forall i \in I \\ & \sum_{j \in J} y_{j} = k & \\ & z_{i} \in \{ 0,1 \} & \forall i \in I \\ & y_j \in \{ 0,1 \} & \forall j \in J \end{array} $$目的関数は，被覆されない顧客の数を最小化することを表す． 最初の制約は，顧客 $i$ が $S$ 内の点にグラフ $G_{\theta}$ 上で隣接するか，そうでない場合には 変数 $z_i$ が $1$ になることを規定する． 2番目の制約は，開設された施設の数が $k$ 個であることを規定する．
この問題を $G_{\theta}$ 上の $k$点被覆問題と呼ぶ． $k$点被覆問題を部分問題として用いて最適な $\theta$ （上の問題の最適値が $0$ になる最小の $\theta$）は2分探索で求めることができる．
def kcover(I, J, c, k):
"""kcover -- minimize the number of uncovered
customers from k facilities.
Parameters:
- I: set of customers
- J: set of potential facilities
- c[i,j]: cost of servicing customer i from facility j
- k: number of facilities to be used
Returns a model, ready to be solved.
"""
model = Model("k-center")
z, y, x = {}, {}, {}
for i in I:
z[i] = model.addVar(vtype="B", name="z(%s)" % i)
for j in J:
y[j] = model.addVar(vtype="B", name="y(%s)" % j)
for i in I:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in I:
model.addConstr(quicksum(x[i, j] for j in J) + z[i] == 1, "Assign(%s)" % i)
for j in J:
model.addConstr(x[i, j] <= y[j], "Strong(%s,%s)" % (i, j))
model.addConstr(quicksum(y[j] for j in J) == k, "k_center")
model.setObjective(
quicksum(z[i] for i in I) + 0.001 * quicksum(c[i, j] * x[i, j] for (i, j) in x),
GRB.MINIMIZE,
)
model.update()
model.__data = x, y, z
return model
def solve_kcenter(I, J, c, k, delta):
"""solve_kcenter -- locate k facilities minimizing
distance of most distant customer.
Parameters:
I - set of customers
J - set of potential facilities
c[i,j] - cost of servicing customer i from facility j
k - number of facilities to be used
delta - tolerance for terminating bisection
Returns:
- list of facilities to be used
- edges linking them to customers
"""
model = kcover(I, J, c, k)
x, y, z = model.__data
facilities, edges = [], []
LB = 0
UB = max(c[i, j] for (i, j) in c)
while UB - LB > delta:
theta = (UB + LB) / 2.0
# print "\n\ncurrent theta:", theta
for j in J:
for i in I:
if c[i, j] > theta:
x[i, j].UB = 0
else:
x[i, j].UB = 1.0
model.update()
model.Params.OutputFlag = 0 # silent mode
model.Params.Cutoff = 0.1
model.optimize()
if model.status == GRB.Status.OPTIMAL:
UB = theta
facilities = [j for j in y if y[j].X > 0.5]
edges = [(i, j) for (i, j) in x if x[i, j].X > 0.5]
else: # infeasibility > 0:
LB = theta
return facilities, edges
random.seed(67)
n = 200
I, J, c, x_pos, y_pos = make_data(n)
k = 5
delta = 1.0e-4
facilities, edges = solve_kcenter(I, J, c, k, delta)
print("Selected facilities:", facilities)
print(
"Max distance from a facility to a customer: ", max([c[i, j] for (i, j) in edges])
)
draw_center(facilities, edges, x_pos, y_pos)
