집합 코팅 문제를 PuLP로 해결

집합 피복 문제를 Python으로 풀고 싶다고 생각했을 때 조사하고 있다고 이 사이트를 만났다.
h tp // w w. 응 ct9. 네. jp / m_ 히로이 / ぃght / 푸 lp06. 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

무사히 실행 결과를 그릴 수 있었다.