集合被覆問題をPuLPで解く

集合被覆問題をPythonで解きたい、と思ったときに調べているとこちらのサイトに出会った。
http://www.nct9.ne.jp/m_hiroi/light/pulp06.html
このサイトにある集合被覆問題についてのコードは私が調査していた問題にちょうどマッチしていたので参考にさせてもらった。こちらのコードではtkinterで実行結果を描画しているが、私はmatplotlibで表現したいと思い、すこし手を加えさせてもらった。

施設配置問題における集合被覆問題は、施設設置可能集合Mと、利用者集合Nが与えられた時、各施設がなるべくすべての利用者をカバーするように、コストと施設数を最小化し、施設の配置を決定することを目的とする。
定式化すると、

Min\quad \sum_{j=1}^{n}x_{j}\\
s.t.\quad \sum_{j=1}^{n}w_{ij}x_{j}\geqq1\quad \forall i\in M \\
x_{j}, w_{ij}\in 0, 1 \quad \forall j \in N\\
M=\{1,2,3,\cdots,m\}\\
N=\{1,2,3,\cdots,n\}

i : 利用者のいる場所
j：施設設置可能地点
M：利用者集合
N：施設設置可能地点集合
xi：地点 i に設置するか否か
wij：地点 i にいる利用者が施設設置可能地点 j の施設への需要があるか否か

今回は、利用者を平面上にランダムに n 個配置し、半径 r 内にある施設を利用できる条件で、集合被覆問題を解く。
n=25, r=20とし、フィールドを80x80の中で最小の配置を考えることにする。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import random, math, time
import pulp

# 入力データ作成、利用者・施設設置可能点生成
def make_data(n, upper):
    random.seed(5)
    X = [random.randint(0, upper) for i in range(n)]
    Y = [random.randint(0, upper) for i in range(n)]
    return [X, Y]

# 利用者(p1)と施設(p2)の距離を測る
def distance(p1, p2):
    dx = p1[0] - p2[0]
    dy = p1[1] - p2[1]
    return math.sqrt(dx * dx + dy * dy)

# 需要変数(二値変数)行列を生成
def make_matrix(ps, r):
    def check(a):
        if a <= r: # 半径がr以下であれば要素を1に、そうでなければ要素を0
            return 1
        else:
            return 0

    unit_list = [] # plotできるようにn次正方行列に整形
    for n,m in zip(ps[0], ps[1]):
        unit = (n,m)
        unit_list.append(unit)

    return [[check(distance(p1, p2)) for p2 in unit_list] for p1 in unit_list]

# グラフ上に描画する
def draw(ps, fs, r):
    fig = plt.figure(figsize=(7,7))
    ax = plt.axes()
    x = np.arange(0, 80, 0.1)
    plt.scatter(ps[0], ps[1], s=20) # 利用者

    for k, (i,j) in zip(fs, zip(ps[0], ps[1])):
        if k.value() == 1: # 施設を設置する場合
            plt.scatter(i, j, s=20, color='red')
            c = patches.Circle(xy=(i, j), radius=r, fill=False) # カバー範囲
            ax.add_patch(c)
        else:
            pass

    plt.axis('scaled')
    ax.set_aspect('equal') # 縦横比を調整
    plt.grid()
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_xticks(np.linspace(0, 80, 9))
    ax.set_yticks(np.linspace(0, 80, 9))
    plt.show()

# 解く
def solver_set(ps, r):
    size = len(ps[0])
    cs = make_matrix(ps, r) # 入力データ生成
    prob = pulp.LpProblem('set-cover')
    fs = [pulp.LpVariable('f{}'.format(i), cat='Binary') for i in range(size)] # 施設設置変数
    prob += pulp.lpSum(fs) # 目的関数
    for c in cs:
        prob += pulp.lpDot(c, fs) >= 1 # 制約条件

    s = time.time() # 時間計測
    status = prob.solve() # 問題が解けたのか、状態確認
    print('Status', pulp.LpStatus[status])
    print('z 施設数：', prob.objective.value())
    print('処理時間：', time.time() - s)
    draw(ps, fs, r)

n = 25 # 利用者・施設を合わせた地点数
r = 20 # 半径
upper = 80 # 乱数の上限値
solver_set(make_data(n, upper), r)

# 結果
# Status Optimal
# z 施設数： 9.0
# 処理時間： 0.05631732940673828

無事、実行結果を描画できた。