Prepare what you need.

python

import random, numpy as np, pandas as pd, networkx as nx, matplotlib.pyplot as plt
from itertools import chain, combinations
from pulp import *

def draw(g):
    """drawing"""
    nx.draw_networkx_labels(g, pos=pos)
    nx.draw_networkx_nodes(g, node_color='w', pos=pos)
    nx.draw_networkx_edges(g, pos=pos)
    plt.show()
def addvar(cnt=[0], *args, **kwargs):
    """Variable creation"""
    cnt[0] += 1
    return LpVariable('v%d'%cnt[0], lowBound=0, *args, **kwargs)

Let's create a problem at random. Demand was set only between some nodes.

python

n = 16 #Number of nodes
g = nx.random_graphs.fast_gnp_random_graph(n, 0.26, 8) #Graph
rn = g.nodes() #Node list
pos = nx.spring_layout(g, pos={i:(i/4, i%4) for i in rn}) #Node position
for i, j in g.edges():
    v = pos[i] - pos[j]
    g[i][j]['dist'] = np.sqrt(v.dot(v)) #distance
dems = random.sample(list(combinations(rn, 2)), 10) #demand
draw(g)

I will try to find the minimum spanning tree.

python

h = nx.minimum_spanning_tree(g, 'dist')
draw(h)

I will try to solve it.

python

rved = [(j, i) for i, j in g.edges()] #Reverse side
#Flow rate per demand
a = pd.DataFrame([(df, dt, ef, et, g.edge[ef][et]['dist'], addvar()) 
                  for df, dt in dems for ef, et in chain(g.edges(), rved)], 
                 columns=['DeFr', 'DeTo', 'EdFr', 'EdTo', 'Dist', 'Var'])
#Whether to install the side
b = pd.DataFrame([(fr, to, g.edge[fr][to]['dist'], addvar(cat=LpBinary)) 
                  for fr, to in g.edges()], columns=['Fr', 'To', 'Dist', 'Var'])
res = [] #Solution(setup cost,Convenience,Graph)
mxcst = 999 #Installation cost upper limit
while True:
    m = LpProblem() #Mathematical model
    m += lpDot(a.Dist, a.Var) + lpDot(b.Dist, b.Var)*1e-6 #Objective function(Convenience)
    m += lpDot(b.Dist, b.Var) <= mxcst #Installation cost upper limit
    for _, r in a.iterrows():
        i, j = r.EdFr, r.EdTo
        if i > j: i, j = j, i
        #Install when flowing the flow rate
        m += r.Var <= lpSum(b.query('Fr==%s & To==%s'%(i,j)).Var)
    for (df, dt), c in a.groupby(('DeFr', 'DeTo')):
        for nd in rn:
            z = 1 if nd == df else -1 if nd == dt else 0
            #Flow rate storage
            m += lpSum(c.query('EdFr == %s'%nd).Var) == \
                 lpSum(c.query('EdTo == %s'%nd).Var) + z
    m.solve() #Solution
    if m.status != 1: break
    a['Val'] = a.Var.apply(lambda v: value(v)) #result(Flow rate)
    b['Val'] = b.Var.apply(lambda v: value(v)) #result(Whether to install)
    cst = value(lpDot(b.Dist, b.Var)) #setup cost
    val = value(m.objective) #Convenience
    mxcst = cst - 1e-3 #Next installation cost upper limit
    h = nx.Graph()
    h.add_edges_from([(r.Fr, r.To) for _, r in b[b.Val == 1].iterrows()])
    res.append((cst, val, h))

Let's see the result. The higher the installation cost, the better the convenience.

python

plt.rcParams['font.family'] = 'IPAexGothic'
plt.plot([r[0] for r in res], [r[1] for r in res])
plt.xlabel('setup cost')
plt.ylabel('Convenience')

It is a graph when the convenience is the best.

python

draw(res[0][2])

I tried to make a GIF animation of the transition of installation cost.

that's all