施設配置問題としての被覆問題と Gurobi による実装 ref: https://qiita.com/osaphex/items/b420f72130fffb58b4e2

file0.tex

\text{Min} \quad \sum_{j \in V} x_j \\
\text{s.t.} \quad \sum_{j \in V_i} x_j \geq 1 \qquad \forall i \in V \\
\qquad x_j \in 0,1 \qquad \forall j \in V

file1.tex

\text{Max} \quad \sum_{i \in V} a_i y_i \\
\text{s.t.} \quad \sum_{j \in V_i} x_j \geq y_i \qquad \forall i \in V \\
\qquad \sum_{j \in V} x_j = p \\
\qquad x_j \in 0, 1 \qquad \forall j \in V \\
\qquad y_i \in 0, 1 \qquad \forall i \in V

file2.py

client = {
    0: {"name": "Client 0", "covered by": [0, 1, 4]},
    1: {"name": "Client 1", "covered by": [0]},
    2: {"name": "Client 2", "covered by": [2, 3, 4]},
    3: {"name": "Client 3", "covered by": [2, 5]},
    4: {"name": "Client 4", "covered by": [2, 5]},
    5: {"name": "Client 5", "covered by": [0, 3]},
    6: {"name": "Client 6", "covered by": [2, 3, 4]},
    7: {"name": "Client 7", "covered by": [1, 4]},
    8: {"name": "Client 8", "covered by": [1, 4, 5]},
}

facility = {
    0: {"name": "Facility 0", "cover": [0, 1, 5]},
    1: {"name": "Facility 1", "cover": [0, 7, 8]},
    2: {"name": "Facility 2", "cover": [2, 3, 4, 6]},
    3: {"name": "Facility 3", "cover": [2, 5, 6]},
    4: {"name": "Facility 4", "cover": [0, 2, 6, 7, 8]},
    5: {"name": "Facility 5", "cover": [3, 4, 8]},
}

file3.py

from gurobipy import *

# モデルの構築
m = Model()

# 決定変数
x = {}
for j in range(len(facility)):
    x[j] = m.addVar(vtype=GRB.BINARY, name="x%d" % j)

m.update()

# 制約条件
for i in range(len(client)):
    m.addConstr( quicksum(x[j] for j in client[i]["covered by"]) >= 1 )

# 目的関数
m.setObjective( quicksum( x[j] for j in range(len(facility)) ), GRB.MINIMIZE )

# 計算の実行
m.optimize()

file4.py

import matplotlib.pyplot as plt
import networkx as nx

fig = plt.figure(figsize=(6, 6))
ax = plt.subplot(111)
ax.set_title("SCLP result", fontsize=18)

G = nx.DiGraph()

for key_client in client:
    G.add_node(client[key_client]["name"])

for key_facility in facility:
    if x[key_facility].x == 1:
        G.add_node(facility[key_facility]["name"])

# 計算結果となるエッジを追加
for key_facility in facility:
    if x[key_facility].x == 1:
        for i in facility[key_facility]["cover"]:
            G.add_edge(client[i]["name"], facility[key_facility]["name"])

nx.draw(G, with_labels=True)
plt.savefig("sclp_result.png", format="PNG")
plt.show()

file5.py

client = {
    0: {"name": "Client 0", "covered by": [0, 1, 4]},
    1: {"name": "Client 1", "covered by": [0]},
    2: {"name": "Client 2", "covered by": [2, 3, 4]},
    3: {"name": "Client 3", "covered by": [2, 5]},
    4: {"name": "Client 4", "covered by": [2, 5]},
    5: {"name": "Client 5", "covered by": [0, 3]},
    6: {"name": "Client 6", "covered by": [2, 3, 4]},
    7: {"name": "Client 7", "covered by": [1, 4]},
    8: {"name": "Client 8", "covered by": [1, 4, 5]},
}

facility = {
    0: {"name": "Facility 0", "cover": [0, 1, 5]},
    1: {"name": "Facility 1", "cover": [0, 7, 8]},
    2: {"name": "Facility 2", "cover": [2, 3, 4, 6]},
    3: {"name": "Facility 3", "cover": [2, 5, 6]},
    4: {"name": "Facility 4", "cover": [0, 2, 6, 7, 8]},
    5: {"name": "Facility 5", "cover": [3, 4, 8]},
}

p = 2

file6.py

from gurobipy import *

# モデルの構築
m = Model()

# 決定変数
x = {}
y = {}

for j in range(len(facility)):
    x[j] = m.addVar(vtype=GRB.BINARY, name="x%d" % j)

for i in range(len(client)):
    y[i] = m.addVar(vtype=GRB.BINARY, name="y%d" % i)

m.update()

# 制約条件
for i in range(len(client)):
    m.addConstr( quicksum(x[j] for j in client[i]["covered by"]) >= y[i] )

m.addConstr( quicksum(x[j] for j in range(len(facility)) ) == p )

# 目的関数
m.setObjective( quicksum(client[i]["population"] * y[i] for i in range(len(client))) , GRB.MAXIMIZE)

# 計算の実行
m.optimize()

file7.py

import matplotlib.pyplot as plt
import networkx as nx

fig = plt.figure(figsize=(6, 6))
ax = plt.subplot(111)
ax.set_title("MCLP result", fontsize=18)

G = nx.DiGraph()

for key_client in client:
    G.add_node(client[key_client]["name"])

for key_facility in facility:
    if x[key_facility].x == 1:
        G.add_node(facility[key_facility]["name"])

for key_facility in facility:
    if x[key_facility].x == 1:
        for i in facility[key_facility]["cover"]:
            if y[i].x == 1:
                G.add_edge(client[i]["name"], facility[key_facility]["name"])

nx.draw(G, with_labels=True)
plt.savefig("mclp_result.png", format="PNG")
plt.show()