I’ve begun learning Python and Gurobi. However, I seem to encounter an issue with my code. Instead of getting the expected solution of 9, I’m getting 49, as the code sums up over every arc. I’m having trouble pinpointing where exactly the problem lies.

import gurobipy as gp
from gurobipy import GRB, Model

import matplotlib.pyplot as plt
import networkx as nx

#item sizes
a = [2, 1, 3]

# profits
p = [5, 4, 3]

# knapsack capacity
b = 5

n = len(a)



model=Model("knapsackL")
model.ModelSense =GRB.MAXIMIZE

weights = {}
#Assign weights for vertical arcs
for i in range(1, b+1):
    for j in range(n+1):
        weights[(i-1, j), (i, j)] = 0

# Assign weights for horizental arcs
for i in range(b+1):
    for j in range(1,n+1):
        weights[(i, j-1), (i, j)] = 0

# Assign weights diagonal arcs
for j in range(1, n+1):
    for i in range(b - a[j-1] + 1):
        weights[(i, j - 1), (i + a[j-1], j)] = p[j-1]

# Filter the outgoing and the ingoing arcs on each node
outgoing = []
ingoing = []
for i in range(b+1):
    for j in range(n+1):
        outgoing += [arc for arc in weights.keys() if (i, j) == arc[0] and not ((0, 0) in arc or (b, n) in arc)]
        ingoing += [arc for arc in weights.keys() if (i, j) == arc[1] and not ((0, 0) in arc or (b, n) in arc)]

# Create binary decision variables for each arc
x={}
for arc, profit in weights.items():
    x[arc]=model.addVar(vtype='b', obj= profit)

#Constraint to ensure exactly one unit of flow leaves node (0,0)
model.addConstr(x[((0,0), (1,0))] + x[((0,0), (0,1))] + x[((0,0), (a[0],1))] == 1)

#Constraint to ensure exactly one unit of flow enters node (b,n)
model.addConstr(x[((b,n-1), (b,n))] + x[((b-1,n), (b,n))] + x[((b-a[n-1],n-1), (b,n))] == 1)

#Flow balance constraint for every other node
model.addConstr(gp.quicksum(x[arc] for arc in outgoing) == gp.quicksum(x[arc] for arc in ingoing))

model.optimize()
#Create a graph object for a visualisation help
G = nx.Graph()

#Add nodes and edges to the graph
for edge, weight in weights.items():
    G.add_edge(edge[0], edge[1], weight=weight)

# Plot the graph
pos = {(x, y): (y, -x) for x, y in G.nodes()}  # positions for all nodes

# Draw nodes
nx.draw_networkx_nodes(G, pos, node_size=700)

# Draw edges
nx.draw_networkx_edges(G, pos)

# Draw edge labels with adjusted positions
edge_labels = {(u, v): f"{d['weight']}" for u, v, d in G.edges(data=True)}

nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, label_pos=0.4)
# Draw node labels

nx.draw_networkx_labels(G, pos)

# Get the current figure
fig = plt.gcf()

# Extract selected edges
selected_edges = [edge for edge in x.keys() if x[edge].x > 0.5]  # Assuming x[edge] > 0.5 means the edge is selected

# Draw selected edges with a different color
nx.draw_networkx_edges(G, pos, edgelist=selected_edges, edge_color='red')

plt.title("Graph in 2D grid")

plt.show()

when I print outgoing, ingoing, and x, I get

outgoing = [((0, 1), (1, 1)), ((0, 1), (0, 2)),…….]

ingoing = [((0, 1), (0, 2)), ((0, 2), (0, 3)), ((0, 1), (1, 1)), ((1, 0), (1, 1)), ((0, 2), (1, 2)), ((1, 1), (1, 2)), ((0, 1), (1, 2)), ((0, 3), (1, 3)), ((1, 2), (1, 3)), ((1, 0), …….]

x = {((0, 0), (1, 0)): <gurobi.Var Awaiting Model Update>, …………}

so the constraint

model.addConstr(gp.quicksum(x[arc] for arc in outgoing) == gp.quicksum(x[arc] for arc in ingoing))

should work !!!

the outgoing list contains the leaving arcs except (0, 0) and (b, n)
the ingoing list contains the entring arcs except (0, 0) in arc or (b, n)
then this constraint iterate on both lists

model.addConstr(gp.quicksum(x[arc] for arc in outgoing) == gp.quicksum(x[arc] for arc in ingoing))