施設配置問題 (Facility Location Problem) と Gurobi による実装 ref: https://qiita.com/osaphex/items/ea20f17855ca25da402e

file0.tex

\text{Max} \quad \sum_i \sum_j c_{ij} x_{ij} + \sum_i f_i y_i \\
\text{s.t.} \quad \sum_i x_{ij} = 1 \quad \forall j \\
\qquad x_{ij} \leq y_i \quad \forall i,j \\
\qquad x_{ij} \geq 0 \quad \forall i,j \\
\qquad y_i \in {0, 1} \quad \forall i,j

file1.tex

\text{Min} \quad \sum_{i,j} h_i d_{ij} x_{ij}\\
\text{s.t.} \quad \sum_j x_{ij} = 1 \qquad \forall i \\
\qquad x_{ij} \leq y_{j} \qquad \forall i, j \\
\qquad \sum_j y_j = p \\
\qquad x_{ij}, y_j \in 0,1 \qquad \forall i, j

file10.py

from gurobipy import *
import math

# 距離
def distance(a,b):
    dx = a[0] - b[0]
    dy = a[1] - b[1]
    return math.sqrt(dx*dx + dy*dy)

d = {}
for i in range(len(client)):
    for j in range(len(facility)):
        d[(i,j)] = distance(client[i]["pos"], facility[j]["pos"])

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

# 決定変数
W = m.addVar(lb=0, vtype=GRB.CONTINUOUS, name="W")

x = {}
y = {}

for i in range(len(client)):
    for j in range(len(facility)):
        x[(i,j)] = m.addVar(lb=0, vtype=GRB.CONTINUOUS, name="x%d,%d" % (i,j))

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

m.update()

# 制約条件
for i in range(len(client)):
    for j in range(len(facility)):
        m.addConstr(x[(i,j)] <= y[j])

for i in range(len(client)):
    m.addConstr(quicksum(x[(i,j)] for j in range(len(facility))) == 1)

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

for i in range(len(client)):
    m.addConstr(W >= quicksum(client[i]["demand"] * d[(i, j)] * x[(i, j)] for j in range(len(facility)) ))

# 目的関数
m.setObjective(W)

# 計算の実行
m.optimize()

file11.py

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

G = nx.DiGraph()

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

for key_facility in facility:
    G.add_node(facility[key_facility]["name"])

pos = {}
for key_client in client:
    pos[client[key_client]["name"]] = client[key_client]["pos"]

for key_facility in facility:
    pos[facility[key_facility]["name"]] = facility[key_facility]["pos"]

# 計算結果のエッジを追加
for i in range(len(client)):
    for j in range(len(facility)):
        if x[(i, j)].x == 1:
            G.add_edge(client[i]["name"], facility[j]["name"])

nx.draw(G, pos, with_labels=True)
plt.savefig("p-center_result.png", format="PNG")
plt.show()

file2.tex

\text{Min} \quad W \\
\text{s.t.} \quad \sum_{j} x_{ij} = 1 \qquad \forall i \\
\qquad \sum_{j} y_j = p \\
\qquad x_{ij} \leq y_j \qquad \forall i,j \\
\qquad W \geq \sum_{j} h_i d_{ij} x_{ij} \qquad \forall i \\
\qquad y_j \in 0, 1 \qquad \forall j \\
\qquad x_{ij} \geq 0 \qquad \forall i,j

file3.py

client = {
    0: {"name": "Client 0", "pos": (0, 1.5), "demand": 1},
    1: {"name": "Client 1", "pos": (2.5, 1.2), "demand": 1}
}

facility = {
    0: {"name": "Facility 0", "pos" : (0, 0), "charge" : 3},
    1: {"name": "Facility 1", "pos" : (0, 1), "charge" : 2},
    2: {"name": "Facility 2", "pos" : (0, 2), "charge" : 3},
    3: {"name": "Facility 3", "pos" : (1, 0), "charge" : 2},
    4: {"name": "Facility 4", "pos" : (1, 1), "charge" : 3},
    5: {"name": "Facility 5", "pos" : (1, 2), "charge" : 3},
    6: {"name": "Faciliyt 6", "pos" : (2, 0), "charge" : 4},
    7: {"name": "Facility 7", "pos" : (2, 1), "charge" : 3},
    8: {"name": "Facility 8", "pos" : (2, 2), "charge" : 2},
}

file4.py

import matplotlib.pyplot as plt
import networkx as nx

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

G = nx.DiGraph()

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

for key_facility in facility:
    G.add_node(facility[key_facility]["name"])

pos = {}
for key_client in client:
    pos[client[key_client]["name"]] = client[key_client]["pos"]

for key_facility in facility:
    pos[facility[key_facility]["name"]] = facility[key_facility]["pos"]

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

file5.py

from gurobipy import *
import math

# 距離
def distance(a,b):
    dx = a[0] - b[0]
    dy = a[1] - b[1]
    return math.sqrt(dx*dx + dy*dy)

d = {}
for i in range(len(client)):
    for j in range(len(facility)):
        d[(i,j)] = distance(client[i]["pos"], facility[j]["pos"])

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

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

for i in range(len(client)):
    for j in range(len(facility)):
        x[(i,j)] = m.addVar(lb=0, vtype=GRB.CONTINUOUS, name="x%d,%d" % (i,j))

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

m.update()

# 制約条件
for i in range(len(client)):
    for j in range(len(facility)):
        m.addConstr(x[(i,j)] <= y[j])

for i in range(len(client)):
    m.addConstr(quicksum(x[(i,j)] for j in range(len(facility))) == 1)

# 目的関数
m.setObjective( quicksum(facility[j]["charge"]*y[j] + quicksum(client[i]["demand"]*d[(i,j)]*x[(i,j)]
                for i in range(len(client))) for j in range(len(facility))) )

# 計算の実行
m.optimize()

file6.py

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

G = nx.DiGraph()

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

for key_facility in facility:
    G.add_node(facility[key_facility]["name"])

pos = {}
for key_client in client:
    pos[client[key_client]["name"]] = client[key_client]["pos"]

for key_facility in facility:
    pos[facility[key_facility]["name"]] = facility[key_facility]["pos"]

# 計算結果のエッジを追加
for i in range(len(client)):
    for j in range(len(facility)):
        if x[(i, j)].x == 1:
            G.add_edge(client[i]["name"], facility[j]["name"])

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

file7.py

client = {
    0: {"name": "Client 0", "pos": (0, 1.5), "demand": 1},
    1: {"name": "Client 1", "pos": (2.5, 1.2), "demand": 1}
}

facility = {
    0: {"name": "Facility 0", "pos" : (0, 0), "charge" : 3},
    1: {"name": "Facility 1", "pos" : (0, 1), "charge" : 2},
    2: {"name": "Facility 2", "pos" : (0, 2), "charge" : 3},
    3: {"name": "Facility 3", "pos" : (1, 0), "charge" : 2},
    4: {"name": "Facility 4", "pos" : (1, 1), "charge" : 3},
    5: {"name": "Facility 5", "pos" : (1, 2), "charge" : 3},
    6: {"name": "Faciliyt 6", "pos" : (2, 0), "charge" : 4},
    7: {"name": "Facility 7", "pos" : (2, 1), "charge" : 3},
    8: {"name": "Facility 8", "pos" : (2, 2), "charge" : 2},
}

p = 2

file8.py

from gurobipy import *
import math

# 距離
def distance(a,b):
    dx = a[0] - b[0]
    dy = a[1] - b[1]
    return math.sqrt(dx*dx + dy*dy)

d = {}
for i in range(len(client)):
    for j in range(len(facility)):
        d[(i,j)] = distance(client[i]["pos"], facility[j]["pos"])

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

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

for i in range(len(client)):
    for j in range(len(facility)):
        x[(i,j)] = m.addVar(lb=0, vtype=GRB.CONTINUOUS, name="x%d,%d" % (i,j))

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

m.update()

# 制約条件
for i in range(len(client)):
    for j in range(len(facility)):
        m.addConstr(x[(i,j)] <= y[j])

for i in range(len(client)):
    m.addConstr(quicksum(x[(i,j)] for j in range(len(facility))) == 1)

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

# 目的関数
m.setObjective( 
    quicksum( quicksum(client[i]["demand"] * d[(i, j)] * x[(i, j)] 
                       for i in range(len(client)) ) for j in range(len(facility)) )
)

# 計算の実行
m.optimize()

file9.py

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

G = nx.DiGraph()

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

for key_facility in facility:
    G.add_node(facility[key_facility]["name"])

pos = {}
for key_client in client:
    pos[client[key_client]["name"]] = client[key_client]["pos"]

for key_facility in facility:
    pos[facility[key_facility]["name"]] = facility[key_facility]["pos"]

# 計算結果のエッジを追加
for i in range(len(client)):
    for j in range(len(facility)):
        if x[(i, j)].x == 1:
            G.add_edge(client[i]["name"], facility[j]["name"])

nx.draw(G, pos, with_labels=True)
plt.savefig("p-median_result.png", format="PNG")
plt.show()