Use Gurobi to create networkx.Graph that has highest edge connectivity

Question

I have the following graph G

G is created using the following code

import networkx as nx
import matplotlib.pyplot as plt

G = nx.hoffman_singleton_graph()
pos = nx.spring_layout(G)
nx.draw(G, pos=pos)
nx.draw_networkx_labels(G=G, pos=pos)
plt.show()

G consists of 50 nodes. I only want to include 25 nodes. In addition, I only want to include the nodes (and edges) that maximizes the connectivity between node A (=5) and node B (=20).

I wrote the following code:

import numpy as np
import gurobipy as grb
from networkx.algorithms.connectivity import local_edge_connectivity


A = 5
B = 20
nodes = list(G.nodes)
n_nodes = len(nodes)
edges = nx.to_numpy_array(G, nodelist=nodes)
thresh_nodes = 25


model = grb.Model()
F = model.addMVar(n_nodes, vtype=grb.GRB.BINARY)

model.addConstr(F.sum() == thresh_nodes)
model.addConstr(F[nodes.index(A)] == 1)
model.addConstr(F[nodes.index(B)] == 1)


E = F * edges

model.setObjective(local_edge_connectivity(nx.from_numpy_array(A=E), A, B), grb.GRB.MAXIMIZE)
model.optimize()

This results in an error because nx.from_numpy_array() cannot deal the datatype of E. How do I create a temporary np.ndarray of the .X-values ({0, 1}) of E and use this to determine the solution(s)?

Answer 1

This can be solved as a MIP (Mixed Integer Program) as you are attempting to do. The below is an implementation in pyomo, and I'm sure it can easily be translated direct to gurobipy if desired. It solves very quickly.

Notes:

• local edge connectivity is synonymous with max unit flow, so we can structure this as a max flow with unit capacity
• we need to augment the non-directed graph to make it directed (essentially doubling the edges)
• set up flow balance and flow limits...
• The local connectivity is 7, which can also be achieved by chopping the threshold down to 15

Code:

import pyomo.environ as pyo
import networkx as nx
import matplotlib.pyplot as plt

G = nx.hoffman_singleton_graph()

edges = list(G.edges())
# augment the graph with reverse edges
edges.extend([(n1, n0) for (n0, n1) in edges])
node_edges = {n:[e for e in edges if e[0]==n or e[1]==n] for n in G.nodes()}

m = pyo.ConcreteModel('graph model')

# SETS
m.N = pyo.Set(initialize=list(G.nodes()))
m.E = pyo.Set(initialize=edges)

# VARS
m.active_node = pyo.Var(m.N, domain=pyo.Binary, doc='Node Active')
m.max_flow = pyo.Var(doc='Max Flow')
m.flow = pyo.Var(m.E, domain=pyo.NonNegativeReals, doc='flow')

s=5
t=20
threshold=25  # <- this could be cleaned up as we know s, t MUST be active

# OBJ: maximize flow
m.obj = pyo.Objective(expr=m.max_flow, sense=pyo.maximize)

# CONSTRAINTS
# 1.  Flow balance
@m.Constraint(m.N)
def flow_balance(m, n):
    inbounds = [e for e in m.E if e[1]==n]
    outbounds = [e for e in m.E if e[0]==n]
    if n == s:
        return m.max_flow <= sum(m.flow[e] for e in outbounds) - sum(m.flow[e] for e in inbounds)
    elif n == t:
        return m.max_flow <= sum(m.flow[e] for e in inbounds) - sum(m.flow[e] for e in outbounds)
    else:
        return sum(m.flow[e] for e in outbounds) == sum(m.flow[e] for e in inbounds)

# 2.  Limit Flow to selection
@m.Constraint(m.E)
def flow_limit_1(m, *e):
    return m.flow[e] <= m.active_node[e[1]]

@m.Constraint(m.E)
def flow_limit_2(m, *e):
    return m.flow[e] <= m.active_node[e[0]]

# 3.  Limit selections
m.node_limit = pyo.Constraint(expr=sum(m.active_node[n] for n in m.N) <= threshold)

# Check it...
m.pprint()

# solve it...
opt = pyo.SolverFactory('appsi_highs')
res = opt.solve(m, tee=True)

m.active_node.display()

print(f'Conectivity: {pyo.value(m.max_flow)}')

# prune it for display
active_nodes = [n for n in m.N if pyo.value(m.active_node[n])]
active_nodes.extend([s, t])
for node in m.N:
    if node not in active_nodes:
        G.remove_node(node)


pos = nx.spring_layout(G)
pos[s] = [-10,0]
pos[t] = [10,0]
nx.draw(G, pos=pos)
nx.draw_networkx_labels(G=G, pos=pos)
plt.show()

Output (partial)