diff --git a/ILP/butadien/butadiensynthesis.py b/ILP/butadien/butadiensynthesis.py
index a0a1c44..85ee175 100644
--- a/ILP/butadien/butadiensynthesis.py
+++ b/ILP/butadien/butadiensynthesis.py
@@ -89,9 +89,9 @@ HYPERGRAPH3 = {
 
 
 FIXED_FLOWS = {
-    #213: 3,
+    213: 3,
     #2: 3,
-    #4: 1,
+    4: 1,
     #23: 1,
     #41: 1,
     #42: 1,
@@ -102,7 +102,9 @@ def build_model(name, hyperedges, elmax, el1, el2, el3, excluded_support=None):
 
     x = {e_id: model.addVar(vtype=GRB.INTEGER, lb = 0, name = f"x_{e_id}") for e_id in hyperedges}
     b = {e_id: model.addVar(vtype=GRB.BINARY, name = f"b_{e_id}") for e_id in hyperedges}
-
+    n1 = {e_id: model.addVar(vtype=GRB.BINARY, name = f"n1_{e_id}") for e_id in hyperedges}
+    n2 = {e_id: model.addVar(vtype=GRB.BINARY, name = f"n2_{e_id}") for e_id in hyperedges}
+    n3 = {e_id: model.addVar(vtype=GRB.BINARY, name = f"n3_{e_id}") for e_id in hyperedges}
  
 
     vertices = set(v for tails, heads in hyperedges.values() for v in tails + heads)
@@ -129,47 +131,38 @@ def build_model(name, hyperedges, elmax, el1, el2, el3, excluded_support=None):
     for e_id in hyperedges:
         model.addGenConstrIndicator(b[e_id], 0, x[e_id] == 0, name = f"unused_implies_zero_{e_id}")
         model.addConstr(x[e_id] >= b[e_id], name = f"used_implies_positive_flow_{e_id}")
+        #Only if an edge has flow, can it contribute to the likelihood
+        model.addConstr(b[e_id] >= n1[e_id], name = f"only_used_contribute_to_nmr1_{e_id}")
+        model.addConstr(b[e_id] >= n2[e_id], name = f"only_used_contribute_to_nmr2_{e_id}")
+        model.addConstr(b[e_id] >= n3[e_id], name = f"only_used_contribute_to_nmr3_{e_id}")
+        #One reaction can only contribute once
+        model.addConstr(n1[e_id] +  n2[e_id] + n3[e_id] <= 1)
+
+    #Only one edge per nmr can contribute
+    model.addConstr(quicksum(n1[e_id] for e_id in hyperedges) == 1)
+    model.addConstr(quicksum(n2[e_id] for e_id in hyperedges) == 1)
+    model.addConstr(quicksum(n3[e_id] for e_id in hyperedges) == 1)
 
     reaction_path = {}
 
 
     if excluded_support:
         model.addConstr(quicksum(b[e_id] for e_id in excluded_support) <= len(excluded_support) - 1, name = "different_hyperedges",)
-    
-    #Multiplizier den node Wert mit infow + outflow
-    model.ModelSense = GRB.MAXIMIZE
 
-    #Multiply edgelikelihood with the edge use boolean
-    #Adapt to have 
-    '''
-    model.setObjectiveN(
-        quicksum(elmax[e_id] * b[e_id] for e_id in hyperedges),
-        index = 0,
-        priority = 2,
-        name = "maximize_nmr_similarity",
-    )  
-    #Minimize the overall flow
-    model.setObjectiveN(
-        quicksum(-1 * x[e_id] for e_id in hyperedges),
-        index=1,
-        priority= 1,
-        name="minimize_used_hyperedges",
-    )
-    '''
+    model.setObjective(quicksum(1000 * (el1[e_id] * n1[e_id] + el2[e_id] * n2[e_id] + el3[e_id] * n3[e_id]) - (1 / elmax[e_id]) * x[e_id] for e_id in hyperedges if elmax[e_id] != 0),GRB.MAXIMIZE)
+
 
     #model.setObjective(quicksum(1000 * elmax[e_id] * b[e_id] - (1 - elmax[e_id]) * x[e_id] for e_id in hyperedges),GRB.MAXIMIZE)
-    model.setObjective(quicksum(1000 * elmax[e_id] * b[e_id] - (1 / elmax[e_id]) * x[e_id] for e_id in hyperedges if elmax[e_id] != 0),GRB.MAXIMIZE)
+    #model.setObjective(quicksum(1000 * elmax[e_id] * b[e_id] - (1 / elmax[e_id]) * x[e_id] for e_id in hyperedges if elmax[e_id] != 0),GRB.MAXIMIZE)
     #model.setObjective(quicksum(1000 * elmax[e_id] * b[e_id] + np.log(elmax[e_id]) * x[e_id] for e_id in hyperedges if elmax[e_id] != 0),GRB.MAXIMIZE)
     #Excluding creation and destruction only three reactions for three nmr
     #model.addConstr(quicksum(b[e_id] for e_id, (tails, heads) in hyperedges.items() if tails != [] and heads != []) == 3) 
-    #Restrict number of used edges to prevent using all available
-    #model.addConstr(quicksum(x[e_id] for e_id in hyperedges) <= 16)
 
     #2 Butadien create first different molecule and it has to be created first:
     startmolecule = ["Butadien"]
 
     #Including it or not changes first and second solution (sometimes flipped)
-    model.addConstr(quicksum(b[e_id] for e_id, (tails, _) in hyperedges.items() if list(set(tails)) == startmolecule) == 1)
+    #model.addConstr(quicksum(b[e_id] for e_id, (tails, _) in hyperedges.items() if list(set(tails)) == startmolecule) == 1)
     #model.addConstr(b[4] + b[7] == 1)
 
     model.addConstr(quicksum(b[e_id] for e_id, (tails, heads) in hyperedges.items() if tails == [] and heads == startmolecule)== 1)