Question: 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?

Answer: 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