python networkx mixed-integer-programming

Generate all paths that consists of specified number of visits of nodes / edges

In a graph/chain there are 3 different states: ST, GRC_i and GRC_j. The following edges between the states exists:

EDGES = [
    # source, target,  name
    ('ST',    'GRC_i', 'TDL_i'),
    ('ST',    'GRC_j', 'TDL_j'),
    ('GRC_i', 'GRC_j', 'RVL_j'),
    ('GRC_j', 'GRC_i', 'RVL_i'),
    ('GRC_j', 'ST',    'SUL_i'),
    ('GRC_i', 'ST',    'SUL_j'),

]

The values for TDL_i, TDL_i, RVL_i and RVL_j are known. The chain always starts in ST and the final state is always known.

I want to infer SUL_i and SUL_j based on possible paths that satisfy the known information.

For example if we have the following information:

RVL_i = 2
RVL_j = 1
TDL_i = 0
TDL_j = 2

and the final position is GRC_i there are two paths that satisfy this criteria:

ST -> TDL_j -> GRC_j -> RVL_i -> GRC_i -> RVL_j -> GRC_j -> SUL_i -> ST -> TDL_j -> GRC_j -> RVL_i -> GRC_i

ST -> TDL_j -> GRC_j -> SUL_i -> ST -> TDL_j -> GRC_j -> RVL_i -> GRC_i -> RVL_j -> GRC_j -> RVL_i -> GRC_i

Because both paths imply that SUL_i = 1 and SUL_j = 0 we conclude that this is the case.

The following relationships are evident:

The number of visits to ST is equal to SUL_i + SUL_j + 1
The number of visits to GRC_i is equal to TDL_i + RVL_i
The number of visits to GRC_j is equal to TDL_j + RVL_j
The upper-bound of SUL_i is the number of visits to GRC_j
The upper-bound of SUL_j is the number of visits to GRC_i
The maximum total number of steps is 2 * (TDL_i + TDL_j + RVL_i + RVL_i)

I was thinking to solve this as a mixed-integer program.

import networkx as nx
import gurobipy as grb
from gurobipy import GRB
from typing import Literal


def get_SUL(TDL_i: int, TDL_j: int, RVL_i: int, RVL_j: int, final_state: Literal['ST', 'GRC_i', 'GRC_j']):
    G = nx.DiGraph()

    G.add_edges_from([
        ('ST', 'GRC_i'),
        ('ST', 'GRC_j'),
        ('GRC_i', 'GRC_j'),
        ('GRC_j', 'GRC_i'),
        ('GRC_j', 'ST'),
        ('GRC_i', 'ST')
    ])

    n_actions = len(list(G.edges()))
    n_states  = len(list(G.nodes()))

    min_N = TDL_i + TDL_j + RVL_i + RVL_i
    max_N = 2 * (TDL_i + TDL_j + RVL_i + RVL_i)

    for N in range(min_N, max_N + 1):
        m = grb.Model()

        SUL_i = m.addVar(lb=0, ub=TDL_j + RVL_j)
        SUL_j = m.addVar(lb=0, ub=TDL_i + RVL_i)

        # actions
        actions = m.addMVar((n_actions, N), vtype=GRB.BINARY)

        m.addConstr(actions[0,:].sum() == TDL_i)
        m.addConstr(actions[1,:].sum() == TDL_j)
        m.addConstr(actions[2,:].sum() == RVL_i)
        m.addConstr(actions[3,:].sum() == RVL_j)
        m.addConstr(actions[4,:].sum() == SUL_i)
        m.addConstr(actions[5,:].sum() == SUL_j)

        m.addConstrs(actions[:,n].sum() == 1 for n in range(N))

        # states
        states = m.addMVar((n_states, N), vtype=GRB.BINARY)

        m.addConstr(states[0,:].sum() == SUL_i + SUL_j + 1)
        m.addConstr(states[0,:].sum() == TDL_i + RVL_i)
        m.addConstr(states[0,:].sum() == TDL_j + RVL_j)

        m.addConstr(states[0,0] == 1)

        if final_state == 'ST':
            m.addConstr(states[0,-1] == 1)
            m.addConstr(states[1,-1] == 0)
            m.addConstr(states[2,-1] == 0)
        elif final_state == 'GRC_i':
            m.addConstr(states[0,-1] == 0)
            m.addConstr(states[1,-1] == 1)
            m.addConstr(states[2,-1] == 0)
        else:
            m.addConstr(states[0,-1] == 0)
            m.addConstr(states[1,-1] == 0)
            m.addConstr(states[2,-1] == 1)

        m.addConstrs(actions[:,n].sum() == 1 for n in range(N))

        # additional constraints

How do I impose that the action- and states variables are in agreement with each other? For example, the first action can only TDL_i or TDL_j because we start in ST. I can obtain the adjacency matrix using nx.to_numpy_array(G) but how should I incorporate this into the model?

Solution

To make things more readable, I will use the following notations:

S, I and J are the nodes
XY is the number of traversals of edge X -> Y

The unknowns of the problem are IS and JS. They must be non-negative.

Case 1: final state is S

During a path, every node is entered and exited the same number of times. For node I, it means:

IS + IJ = SI + JI

For node J:

JS + JI = SJ + IJ

The equations immediately lead to the solution:

IS = SI + JI - IJ
JS = SJ + IJ - JI

There is a special case to consider: if SI and SJ are 0, we can't leave node S, so only the empty path is possible. If IJ or JI is greater than 0, then no solution is possible.

Case 2: final state is I

For node I, the number of entries is one greater than the number of exits:

IS + IJ + 1 = SI + JI

For node J, they are the same:

JS + JI = SJ + IJ

Which gives:

IS = SI + JI - IJ - 1
JS = SJ + IJ - JI

Case 3: final state is J

Node I is entered and exited the same number of times:

IS + IJ = SI + JI

For node J, the number of entries is one greater than the number of exits:

JS + JI + 1 = SJ + IJ

Which gives:

IS = SI + JI - IJ
JS = SJ + IJ - JI - 1

Code

def solve(SI, SJ, IJ, JI, final_state):
    match final_state:
        case "S":
            if SI == 0 and SJ == 0:
                if IJ == 0 and JI == 0:
                    return (0, 0)
                else:
                    return None
            IS = SI + JI - IJ
            JS = SJ + IJ - JI
        case "I":
            IS = SI + JI - IJ - 1
            JS = SJ + IJ - JI
        case "J":
            IS = SI + JI - IJ
            JS = SJ + IJ - JI - 1

    if IS >= 0 and JS >= 0:
        return (IS, JS)
    else:
        return None

def test(SI, SJ, IJ, JI, final_state):
    res = solve(SI, SJ, IJ, JI, final_state)
    inputs = f"SI={SI}, SJ={SJ}, IJ={IJ}, JI={JI}, final_state={final_state}"
    if res is None:
        print(f"{inputs} => no solution")
    else:
        print(f"{inputs} => IS={res[0]}, JS={res[1]}")

test(SI=0, SJ=2, IJ=1, JI=2, final_state="I")
test(SI=1, SJ=0, IJ=9, JI=9, final_state="I")
test(SI=0, SJ=2, IJ=1, JI=1, final_state="I")
test(SI=0, SJ=2, IJ=1, JI=2, final_state="J")
test(SI=0, SJ=2, IJ=1, JI=1, final_state="J")
test(SI=1, SJ=1, IJ=0, JI=2, final_state="J")
test(SI=2, SJ=0, IJ=2, JI=0, final_state="S")
test(SI=2, SJ=0, IJ=2, JI=1, final_state="S")
test(SI=0, SJ=0, IJ=1, JI=1, final_state="S")

Results

The first case is the one given in the question:

SI=0, SJ=2, IJ=1, JI=2, final_state=I => IS=0, JS=1
SI=1, SJ=0, IJ=9, JI=9, final_state=I => IS=0, JS=0
SI=0, SJ=2, IJ=1, JI=1, final_state=I => no solution
SI=0, SJ=2, IJ=1, JI=2, final_state=J => IS=1, JS=0
SI=0, SJ=2, IJ=1, JI=1, final_state=J => IS=0, JS=1
SI=1, SJ=1, IJ=0, JI=2, final_state=J => no solution
SI=2, SJ=0, IJ=2, JI=0, final_state=S => IS=0, JS=2
SI=2, SJ=0, IJ=2, JI=1, final_state=S => IS=1, JS=1
SI=0, SJ=0, IJ=1, JI=1, final_state=S => no solution