P-median問題とWeber問題のコードについて

Question

前提・実現したいこと

最近プログラミングを始めた超初心者です。
Python3(Anaconda Jupyter Notebook) で数理最適化問題を解いています。
エクセルから座標データを入力し，それをP-median問題を用いて拠点の導出を行い，そこからWeber問題用いて中継となる地点を導出したいと考えています。

現段階としてP-median問題のコードはできているんですが，それ以降のWeber問題が分かりません。また，拠点となった地点に対して色を変えたりしたいのですがどうやればいいのでしょうか?

よろしくお願いします。

該当のソースコード

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 
import networkx as nx 
import pulp
import xlrd
import math
import random


#=====Parameter=====
h_omomi = 1
s_omomi = 1
p_omomi = 1
#===================

#エクセルからデータ取得
book = xlrd.open_workbook("***.xlsx")
sheet = book.sheet_by_name("Sheet")

def get_data(sheet,rowx,colx):
    data = []
    for row in range(rowx, sheet.nrows):  #データ起点から
        value = sheet.cell(row, colx).value
        data.append(value)
    return data
        
x = get_data(sheet, 1, 1)                #データの起点はB1
y = get_data(sheet, 1, 2)                #データの起点はC1
H = get_data(sheet, 1, 3)               #データの起点はD1
S = get_data(sheet, 1, 4)               #データの起点はE1
P = get_data(sheet, 1, 5)               #データの起点はF1
    
#リスト化
x_1 = [int(h) for h in x]
y_1 = [int(i) for i in y]
H_1 = [int(j) for j in Hospital]
S_1 = [int(k) for k in School]
P_1 = [int(l) for l in Population]
    
print(x_1)
print(y_1)
print(H_1)
print(S_1)
print(P_1)

pos = list(zip(x_1,y_1))
print(pos)


# 各駅における需要量の計算
w = []
for i in range(len(H_1)):
    T_W = h_omomi * H_1[i] + s_omomi * S_1[i] + p_omomi * P_1[i]
    w.append(T_W)
print(w)

G = nx.Graph()
for i in range(len(x_1)):
        G.add_node(i, w=w[i], pos=pos[i])

def draw_facilities(G):
    pos = nx.get_node_attributes(G,'pos')
    nx.draw_networkx(G, pos=pos, with_labels=False)
    w = nx.get_node_attributes(G, 'w')
    _ = nx.draw_networkx_labels(G, pos=pos, labels=w, font_color='w')
    return

plt.figure(figsize=(10,10))
draw_facilities(G)
plt.xlabel("x")
plt.ylabel("y")
plt.savefig('plot_1.png')

def solve_p_median_by_pulp(G, p):
    pmedian = pulp.LpProblem(name='5-median', sense=pulp.LpMinimize)
    w = nx.get_node_attributes(G, 'w')
    pos = nx.get_node_attributes(G, 'pos')
    x = {}
    y = {}
    c = {}
    for v in G.nodes():
        x[v] = pulp.LpVariable(name='x{0}'.format(v), cat='Binary')
        for u in G.nodes():
            y[(u, v)] = pulp.LpVariable(name='y{0}_{1}'.format(u, v), cat='Binary')
            c[(u, v)] = math.sqrt(sum([(pos[u][i] - pos[v][i])**2 for i in range(2)]))
    pmedian += pulp.LpAffineExpression([(y[(u, v)], w[u] * c[(u, v)]) for u, v in y])
    pmedian += pulp.LpAffineExpression([(x[v], 1) for v in G.nodes()]) == p
    for u in G.nodes():
        pmedian += pulp.LpAffineExpression([(y[(u, v)], 1) for v in G.nodes()]) == 1
    for u, v in y:
        pmedian += y[(u, v)] - x[v] <= 0
    pmedian.solve()
    for v in x:
        G.node[v]['x'] = x[v].value()
    for u, v in y:
        if y[(u, v)].value() == 1:
            G.add_edge(u, v)
    return G

G = solve_p_median_by_pulp(G,5)
plt.figure(figsize=(10,10))
draw_facilities(G)
   
    
plt.xlabel("x")
plt.ylabel("y")
plt.savefig('***.png')

Answer 1

同じ質問がここにあるので参考になると思います．