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

センター問題

センター問題(center problem)とは, 顧客から最も近い施設への距離の「最大値」を最小にするように,グラフ内の点または枝上,または空間内の任意の点から $k$ 個の施設を選択する問題の総称であり, 施設数 $k$ を指定した問題を$k$-センター問題($k$-center problem) とよぶ.

ここでは,点上に施設を配置する場合を考える.

標準定式化

$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]
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 961 rows, 931 columns and 3630 nonzeros
Model fingerprint: 0x8230ce33
Variable types: 1 continuous, 930 integer (930 binary)
Coefficient statistics:
  Matrix range     [3e-02, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 5e+00]
Presolve time: 0.01s
Presolved: 961 rows, 931 columns, 3630 nonzeros
Variable types: 1 continuous, 930 integer (930 binary)
Found heuristic solution: objective 0.7069622
Found heuristic solution: objective 0.6501211
Found heuristic solution: objective 0.6332376

Root relaxation: objective 1.998137e-01, 507 iterations, 0.01 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    0.19981    0  228    0.63324    0.19981  68.4%     -    0s
H    0     0                       0.5951071    0.19981  66.4%     -    0s
H    0     0                       0.4970149    0.19981  59.8%     -    0s
     0     0    0.22202    0  222    0.49701    0.22202  55.3%     -    0s
     0     0    0.22202    0  228    0.49701    0.22202  55.3%     -    0s
H    0     0                       0.4358186    0.22202  49.1%     -    0s
H    0     0                       0.4298363    0.22202  48.3%     -    0s
H    0     0                       0.4044511    0.22221  45.1%     -    0s
     0     0    0.22952    0  193    0.40445    0.22952  43.3%     -    0s
H    0     0                       0.3950400    0.22952  41.9%     -    0s
     0     0    0.23666    0  169    0.39504    0.23666  40.1%     -    0s
     0     0    0.24335    0  161    0.39504    0.24335  38.4%     -    0s
     0     0    0.24335    0  129    0.39504    0.24335  38.4%     -    0s
H    0     0                       0.3700418    0.24335  34.2%     -    0s
H    0     0                       0.3455534    0.24335  29.6%     -    0s
     0     0    0.26065    0  120    0.34555    0.26065  24.6%     -    0s
     0     0    0.26445    0  112    0.34555    0.26445  23.5%     -    0s
     0     0    0.27611    0  119    0.34555    0.27611  20.1%     -    0s
     0     0    0.29390    0  101    0.34555    0.29390  14.9%     -    0s
     0     0    0.29594    0  119    0.34555    0.29594  14.4%     -    0s
     0     0    0.29740    0  124    0.34555    0.29740  13.9%     -    0s
     0     0    0.29740    0  124    0.34555    0.29740  13.9%     -    0s
H    0     0                       0.3352990    0.29740  11.3%     -    0s
     0     0    0.30794    0   32    0.33530    0.30794  8.16%     -    0s
H    0     0                       0.3336685    0.30794  7.71%     -    0s
     0     0     cutoff    0         0.33367    0.33367  0.00%     -    0s

Cutting planes:
  Gomory: 6
  MIR: 5
  Zero half: 1
  Mod-K: 1
  RLT: 6

Explored 1 nodes (3018 simplex iterations) in 0.27 seconds
Thread count was 16 (of 16 available processors)

Solution count 10: 0.333668 0.335299 0.345553 ... 0.595107

Optimal solution found (tolerance 1.00e-04)
Best objective 3.336684909605e-01, best bound 3.336684909605e-01, gap 0.0000%
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])
)
Selected facilities: [89, 93, 131, 133, 196]
Max distance from a facility to a customer:  0.3061769941895321
draw_center(facilities, edges, x_pos, y_pos)