Find max unit flow in Networkx.Graph using Gurobi

Hi all,

I have a `networkx.Graph`-object called G.

import networkx as nx
import gurobipy as grb
from gurobipy import GRB
import matplotlib.pyplot as plt

G = nx.hoffman_singleton_graph()

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

I wrote the following script:


def update_graph(G: nx.Graph, A: int, B: int, thresh: int) -> nx.Graph:
    nodes   = list(G.nodes)
    Adj     = nx.to_numpy_array(G, nodelist=nodes)
    n_nodes = len(nodes)
    edges   = list(zip(*Adj.nonzero()))

    m        = grb.Model()
    N        = m.addMVar(n_nodes, vtype=GRB.BINARY)
    flow     = m.addVars(edges, vtype=GRB.CONTINUOUS, lb=0)
    max_flow = m.addVar(vtype=GRB.INTEGER, lb=0)
    E        = Adj * N

    for n in nodes:
        inbounds  = [(u, v) for u, v in edges if v == n]
        outbounds = [(u, v) for u, v in edges if u == n]

        if n == A:
            m.addConstr(max_flow <= grb.quicksum(flow[e] for e in outbounds) - grb.quicksum(flow[e] for e in inbounds))

        elif n == B:
            m.addConstr(max_flow <= grb.quicksum(flow[e] for e in inbounds) - grb.quicksum(flow[e] for e in outbounds))

        else:
            m.addConstr(grb.quicksum(flow[e] for e in outbounds) == grb.quicksum(flow[e] for e in inbounds))

    for u, v in edges:
        m.addConstr(flow[u, v] <= N[v])
        m.addConstr(flow[u, v] <= N[u])

    m.addConstr(N[nodes.index(A)] == 1)
    m.addConstr(N[nodes.index(B)] == 1)
    m.addConstr(N.sum() <= thresh)

    m.setObjective(max_flow, grb.GRB.MAXIMIZE)
    m.optimize()

    print(f'Connectivity: {max_flow.Xn}')
    n_solutions = m.SolCount

    max_E = 0
    max_s = 0

    for s in range(n_solutions):
        # select the solution with the most edges 
        m.setParam(GRB.Param.SolutionNumber, s)

        E_s = E.sum(axis=1)[
            [nodes.index(A),
             nodes.index(B)]
        ].getValue().sum()

        if E_s > max_E:
            max_s = s
            max_E = E_s

        pool = [n for n in nodes if N[n].Xn > 0.5]
        print(pool)

    m.setParam(GRB.Param.SolutionNumber, max_s)

    for node in nodes:
        if N[node].Xn == 0:
            G.remove_node(node)

    return G


blue = 5
red  = 20
n    = 25

G = update_graph(G=G, A=blue, B=red, thresh=n)

This script almost works, but printing the pool-variable shows that there is something going wrong.. For the second solution there are only 2 active nodes but the unit flow is 7. That's not possible.

By changing the binary values of N the unit flow between the two specified nodes should change. However, there is currently no relationship defined between N and flow. Therefore, the constraints are static.

If N[0] is equal to 0, the first row and column of E should be equal to 0 because the node and hence its edges have "disappeared".

Should E be linked to the flow-constraints? If so, how can I achieve this?

Find max unit flow in Networkx.Graph using Gurobi

Comments

Didn't find what you were looking for?