Gurobi unable to find the most optimal Link Utilization

Question

from gurobipy import Model, GRB, quicksum, Env
def calcMLU(trafficMatrix, adjNodes, OSPFflowPaths, SDNNodes):
  params = {
  "WLSACCESSID": '',
  "WLSSECRET": '',
  "LICENSEID": 123,
  }

  env = Env(params=params)

  # Create a new model
  m = Model("network_utilization", env=env)

  # Hide Output
  m.setParam('OutputFlag', 0)

  links = generateLinks(adjNodes)
  demands = generateDemands(trafficMatrix)
  SDNNodes = SDNNodes
  nonSDNNodes = [n for n in adjNodes.keys() if n not in SDNNodes]
  nextHop = generateNextHop(OSPFflowPaths)
  nodes = list(adjNodes.keys())

  # Capacity for each link in the network
  capacities = {link:1000000000 for link in links}

  # Variables representing the flow on each link for each demand
  flow_vars = m.addVars([(demand, link) for demand in demands for link in links], vtype=GRB.CONTINUOUS, name="flow")

  # Utilization variable for each link
  util_vars = m.addVars(links, vtype=GRB.CONTINUOUS, name="util")

  # Auxiliary variable to represent the maximum utilization
  max_util = m.addVar(vtype=GRB.CONTINUOUS, name="max_util")

  # Objective: Minimize the maximum utilization
  m.setObjective(max_util, GRB.MINIMIZE)

  for demand in demands:
    for node in SDNNodes:  # SDNNodes
      inflow = quicksum(flow_vars[demand, (i, j)] for i, j in links if j == node and (demand, (i, j)) in flow_vars)
      outflow = quicksum(flow_vars[demand, (i, j)] for i, j in links if i == node and (demand, (i, j)) in flow_vars)
      # Source node
      if node == demand[0]:
        m.addConstr(outflow - inflow == demands[demand], f'flow_conservation_source_{demand}_{node}')
        #m.addConstr(inflow == 0, f'flow_conservation_source2_{demand}_{node}')
      # Destination node
      elif node == demand[1]:
        m.addConstr(inflow - outflow == demands[demand], f'flow_conservation_sink_{demand}_{node}')
        #m.addConstr(outflow == 0, f'flow_conservation_sink2_{demand}_{node}')
      # Intermediate nodes
      else:
        m.addConstr(inflow == outflow, f'flow_conservation_transit_{demand}_{node}')


  for demand in demands:
    for node in nonSDNNodes:  # nonSDNNodes
      inflow = quicksum(flow_vars[demand, (i, j)] for i, j in links if j == node and (demand, (i, j)) in flow_vars)
      # If Destination node, outflow must be 0, inflow must be demands
      if node == demand[1]:
        outflow = quicksum(flow_vars[demand, (i, j)] for i, j in links if i == node and (demand, (i, j)) in flow_vars)
        m.addConstr(inflow == demands[demand], f'flow_conservation_sink_{demand}_{node}')
        m.addConstr(outflow == 0, f'flow_conservation_sink2_{demand}_{node}')
        continue

      dest = demand[1]
      next = nextHop[node][dest]
      outflow = flow_vars[demand, (node, next)]
      notOutflow = quicksum(flow_vars[demand, (i, j)] for i, j in links if i == node and (demand, (i, j)) in flow_vars and j != next)

      # Source node
      if node == demand[0]:
        m.addConstr(outflow == demands[demand], f'flow_conservation_source_{demand}_{node}')
        m.addConstr(notOutflow == 0, f'flow_conservation_not_source_{demand}_{node}')
        m.addConstr(inflow == 0, f'flow_conservation_source2_{demand}_{node}')
      # Intermediate nodes
      else:
        m.addConstr(inflow == outflow, f'flow_conservation_transit_{demand}_{node}')
        m.addConstr(notOutflow == 0, f'flow_conservation_not_transit_{demand}_{node}')


  # Capacity and utilization constraints for each link
  for link in links:
      link_flow = quicksum(flow_vars[demand, link] for demand in demands if (demand, link) in flow_vars)
      m.addConstr(link_flow <= capacities[link], f'capacity_{link}')
      m.addConstr(util_vars[link] == link_flow / capacities[link], f'utilization_{link}')

  # Constraints to ensure max_util is at least as large as any utilization
  for link in links:
      m.addConstr(max_util >= util_vars[link], f'max_util_constraint_{link}')


  # Solve the model
  m.optimize()

  # Check the number of constraints
  num_constraints = m.NumConstrs
  print(f"Number of constraints: {num_constraints}")

  # Output the results
  if m.status == GRB.OPTIMAL:
      print("Optimal link utilization:")
      for link in links:
          print(f"{link}: {util_vars[link].X:.2%}")
      print(f"Maximum utilization among all links: {max_util.X:.2%}")
  else:
      print("No optimal solution found.")

  # Transpose into Python Dictionary
  LU = generateLU(links, util_vars)
  MLU = max_util.X

  return MLU

A beginner coder here. Need help on Gurobi, may I know why this optimization produced worse result when more nodes are migrating to SDN. Is there any problem with my formulation, sincerely appreciate for all helps. The flow conservation rules seems got issues but I have no idea which part of it went wrong. The MLU gradually getting worse despite more nodes migrated to SDN.

A fix for the code especially in the flow constraints