diff --git a/python/conda package/k_hop_dominating_set_gurobi/k_hop_dominating_set_gurobi/k_hop_dom_set.py b/python/conda package/k_hop_dominating_set_gurobi/k_hop_dominating_set_gurobi/k_hop_dom_set.py
index e45ffcb5774486d38ea734b00bdea2705d4da3ec..adeecb79973a51d89dd9b644aa8ed948df03417b 100755
--- a/python/conda package/k_hop_dominating_set_gurobi/k_hop_dominating_set_gurobi/k_hop_dom_set.py	
+++ b/python/conda package/k_hop_dominating_set_gurobi/k_hop_dominating_set_gurobi/k_hop_dom_set.py	
@@ -13,20 +13,30 @@ import sys
 import matplotlib.pyplot as plt
 #import lp_to_nx_graph
 from itertools import combinations
+from enum import Enum
 
+class ConnectivityConstraint(Enum):
+    SEPARATORS = 1
+    MILLER_TUCKER_ZEMLIN = 3
+    MARTIN = 4
     
-class DominatingSet:
+class KHopDominatingSet():
     
-    def __init__(self, G, name = "DS"):
+    def __init__(self, G, k, name = "kDS"):
         self.G = G
         self.m = gp.Model(name)
         self.nodes = self.m.addVars(G.nodes, vtype = GRB.BINARY, name = "nodes")
         self.m.setObjective(gp.quicksum(self.nodes), GRB.MINIMIZE)
-    
-        # For each node in the graph the node itself or at least one of it's neighbors has to be part of the dominating set
-        self.m.addConstrs(((gp.quicksum(self.nodes[n] for n in G.neighbors(v)) + self.nodes[v] )>= 1) for v in G.nodes)
+
+        # For each node in the graph the node itself or at least one of it's k-neighbors has to be part of the dominating set
+        for v in G.nodes: 
+            k_neighborhood = set(G.neighbors(v))
+            for i in range(k-1):
+                k_neighborhood.update([w for n in k_neighborhood for w in G.neighbors(n)])
+            self.m.addConstr((gp.quicksum(self.nodes[n] for n in k_neighborhood) +self.nodes[v]) >= 1 )
 
     def solve(self):
+        # TODO: refactor, maybe use self.nodes directly
         ds = {i for i,x_i in enumerate(self.m.getVars()) if x_i.x == 1}
         return ds
     
@@ -40,13 +50,7 @@ class DominatingSet:
         print(f"duration in seconds: {duration_sec}")
         nx.draw_kamada_kawai(self.G, node_color = color_map, with_labels = True)
         plt.show()
-        
-    
-class KHopDominatingSet(DominatingSet):
-    
-    def __init__(self, G, k, name = "kDS"):
-        super().__init__(KHopDominatingSet.make_closure(G,k), name)
-        self.G = G
+
 
     def make_closure(G, k):
         G_prime = nx.Graph(G)
@@ -61,17 +65,78 @@ class KHopDominatingSet(DominatingSet):
 
 class ConnectedKHopDominatingSet(KHopDominatingSet):
     
-    def __init__(self, G, k, name = "CkDS", exclude = {}):
+    def __init__(self, G, k, name = "CkDS", exclude = {}, constraints = ConnectivityConstraint.SEPARATORS):
         self.G = G
         super().__init__(G, k, name)
         
-        if exclude:
-            self.m.addConstrs(self.nodes[v] <= gp.quicksum(self.nodes[w] for w in G.neighbors(v)) for v in G.nodes if v not in exclude)
-        else:
-            self.m.addConstrs(self.nodes[v] <= gp.quicksum(self.nodes[w] for w in G.neighbors(v)) for v in G.nodes)
+        self.constraints = constraints
         
+        if ConnectivityConstraint.SEPARATORS == constraints:
+            if exclude:
+                self.m.addConstrs(self.nodes[v] <= gp.quicksum(self.nodes[w] for w in G.neighbors(v)) for v in G.nodes if v not in exclude)
+            else:
+                self.m.addConstrs(self.nodes[v] <= gp.quicksum(self.nodes[w] for w in G.neighbors(v)) for v in G.nodes)
+      
+        if ConnectivityConstraint.MILLER_TUCKER_ZEMLIN == constraints:
+            # given an undirected graph generates a graph which is bidirectional
+            self.G_prime = self.G.to_directed()
 
-        
+            n = len(self.G.nodes)
+            self.G_prime.add_nodes_from([n+1, n+2])
+            self.G_prime.add_edge(n+1,n+2)
+            initial_nodes = self.G.nodes
+            
+            self.G_prime.add_edges_from((n+1, i) for i in initial_nodes)
+            self.G_prime.add_edges_from((n+2, i) for i in initial_nodes)
+            self.edges = self.m.addVars(self.G_prime.edges, vtype = GRB.BINARY, name="edges")
+            
+            self.auxiliary = self.m.addVars(self.G_prime.nodes, vtype= GRB.INTEGER, name="auxiliary")
+            
+            self.m.addConstr(gp.quicksum(self.edges[(n+2,i)] for i in self.G.nodes) == 1)
+            self.m.addConstrs(gp.quicksum(self.edges[(i,j)] for i in self.G_prime.nodes if i != j and (i,j) in self.G_prime.edges) == 1 for j in self.G.nodes)
+            
+            self.m.addConstrs(self.edges[(n+1,i)] + self.edges[(i,j)] <= 1 for (i,j) in self.G.to_directed().edges)
+
+            self.m.addConstrs((n+1)*self.edges[(i,j)] + self.auxiliary[i] - self.auxiliary[j]+ (n-1)*self.edges[(j,i)] <= n for (i,j) in self.G.to_directed().edges)
+            
+            A_without_E_and_E_prime = [e for e in self.G_prime.edges if e not in self.G.to_directed().edges]
+            self.m.addConstrs((n+1)*self.edges[(i,j)]+self.auxiliary[i]-self.auxiliary[j] <= n for (i,j) in A_without_E_and_E_prime)
+            
+            self.m.addConstr(self.edges[(n+1,n+2)] == 1)
+            
+            self.m.addConstr(self.auxiliary[n+1] == 0)
+            self.m.addConstrs(self.auxiliary[i] >= 1 for i in self.G_prime.nodes if i != n+1)
+            self.m.addConstrs(self.auxiliary[i] <= n+1 for i in self.G_prime.nodes if i != n+1)
+            
+            self.m.addConstrs(self.nodes[i] == 1-self.edges[(n+1,i)] for i in self.G.nodes)
+            
+        if ConnectivityConstraint.MARTIN == constraints:
+            # large integer
+            M = 99999
+            
+            E = [self.G.edges]
+            
+            self.edges = self.m.addVars(((i,j) for i in self.G.nodes for j in self.G.nodes), vtype = GRB.BINARY, name = "edges")
+            
+            self.z = self.m.addVars( (((i,j) ,k) for i in self.G.nodes for j in self.G.nodes for k in self.G.nodes), vtype = GRB.BINARY, name ="z")
+            
+            self.m.addConstr(gp.quicksum(self.edges) == gp.quicksum(self.nodes) - 1)
+            
+            self.m.addConstrs(self.edges[(i,j)] <= self.nodes[i] for (i,j) in self.G.edges)
+            self.m.addConstrs(self.edges[(i,j)] <= self.nodes[j] for (i,j) in self.G.edges)
+            
+            self.m.addConstrs(self.z[((i,j),k)] <= self.edges[(i,j)] for (i,j) in self.G.edges for k in self.G.nodes)
+            self.m.addConstrs(self.z[((j,i),k)] <= self.nodes[k] for (i,j) in self.G.edges for k in self.G.nodes)
+            
+            self.m.addConstrs(self.edges[(i,j)] - M*(3-self.nodes[i]-self.nodes[j]-self.nodes[k]) <= self.z[((i,j),k)] + self.z[((j,i),k)] for i in self.G.nodes for j in self.G.nodes for k in self.G.nodes)
+            self.m.addConstrs(self.z[((i,j),k)] + self.z[((j,i),k)] <= self.edges[(i,j)] +M*(3-self.nodes[i]-self.nodes[j]-self.nodes[k]) for i in self.G.nodes for j in self.G.nodes for k in self.G.nodes)
+                                    
+            self.m.addConstrs(1-M*(2-self.nodes[i]-self.nodes[j]) <= gp.quicksum(self.z[((i,k),j)] for k in self.G.nodes if k != i and k != j if (i,k) in self.G.to_directed().edges)  +self.edges[(i,j)] for i in self.G.nodes for j in self.G.nodes)
+            
+            self.m.addConstrs(self.edges[(i,j)] == 0 for i in self.G.nodes for j in self.G.nodes if (i,j) not in E)
+            self.m.addConstrs(self.z[((i,j),k)] == 0 for i in self.G.nodes for j in self.G.nodes for k in self.G.nodes if (i,j) not in E)
+            
+            
     def min_ij_separator(G, i, j, C_i):
         N_ci = {v for c in C_i for v in G.neighbors(c)}
         G_prime = nx.Graph(G)
@@ -83,12 +148,19 @@ class ConnectedKHopDominatingSet(KHopDominatingSet):
         return R_j.intersection(N_ci)
     
     def solve(self):
-        self.m._vars = self.nodes
-        self.m._G = self.G
-        self.m.Params.lazyConstraints = 1
-        self.m.optimize(RootedConnectecKHopDominatingSet.elim_unconnectivity)    
-        # elif self.constraints == ConnectivityConstraint.INDEGREE:
-        #     self.m.optimize()
+
+        if ConnectivityConstraint.SEPARATORS == self.constraints:
+            self.m._vars = self.nodes
+            self.m._G = self.G
+            self.m.Params.lazyConstraints = 1
+            self.m.optimize(ConnectedKHopDominatingSet.elim_unconnectivity)    
+
+        elif ConnectivityConstraint.INDEGREE == self.constraints:
+            self.m.optimize()
+        elif ConnectivityConstraint.MILLER_TUCKER_ZEMLIN == self.constraints:
+            self.m.optimize()
+        elif ConnectivityConstraint.MARTIN == self.constraints:
+            self.m.optimize()
         
         ds = {i for i,x_i in enumerate(self.m.getVars()) if x_i.x > 0.5}
         return ds
@@ -113,12 +185,16 @@ class ConnectedKHopDominatingSet(KHopDominatingSet):
         
 class RootedConnectecKHopDominatingSet(ConnectedKHopDominatingSet):
     
-    def __init__(self, G, k, root = 0, name = "RCkDS"):
-        super().__init__(G, k, name, exclude = {root})
+    def __init__(self, G, k, root = 0, name = "RCkDS", constraints = ConnectivityConstraint.SEPARATORS):
+        super().__init__(G, k, name, exclude = {root}, constraints = constraints)
         self.root = root
         
         self.m.addConstr(self.nodes[root] >= 1)
         
+        # self.add_all_combinations_root_separators()
+        self.add_all_neighborhood_root_separators()
+        # self.add_single_root_separators()
+        
     def add_single_root_separators(self):
         for i in self.G.nodes:
             min_ij_sep = ConnectedKHopDominatingSet.min_ij_separator(G, i, self.root, {i})
@@ -127,18 +203,39 @@ class RootedConnectecKHopDominatingSet(ConnectedKHopDominatingSet):
         
     # ToDo: refactor!
     def add_all_neighborhood_root_separators(self):
+        number_of_found_separators = 0
         for v in self.G.nodes:
+            # V1
+            to_add = set()
             if v is not self.root and self.root not in self.G.neighbors(v) and v not in self.G.neighbors(self.root) and not set(self.G.neighbors(self.root)).intersection(set(self.G.neighbors(v))):
-                for i in range(2,self.G.degree[v]):
+                for i in range(2,self.G.degree[v]+1):
                     V = {w for w in self.G.neighbors(v)}
                     V.update([v])
                     for i_neighborhood in combinations(V, i):
+                        # V3
+                        # to_add = set()
                         if v in i_neighborhood:
+                            number_of_found_separators += 1
                             min_ij_sep = ConnectedKHopDominatingSet.min_ij_separator(self.G, v, self.root, set(i_neighborhood))
-                            self.m.addConstr(gp.quicksum(self.nodes[s] for s in min_ij_sep) >= self.nodes[v])
+                            # V1
+                            to_add.add(frozenset(min_ij_sep))
+                            
+                            # V2
+                            # self.m.addConstr(gp.quicksum(self.nodes[s] for s in min_ij_sep) >= self.nodes[v])
+                        
+                        # V3
+                        # self.m.addConstrs(gp.quicksum(self.nodes[s] for s in min_ij_sep_to_add) >= self.nodes[w] for w in i_neighborhood for min_ij_sep_to_add in to_add)
+            # V1
+            self.m.addConstrs(gp.quicksum(self.nodes[s] for s in min_ij_sep) >= self.nodes[v] for min_ij_sep in to_add)
+
+ 
     
     def add_all_combinations_root_separators(self):
+        number_of_found_separators = 0
+        cummulated_to_add = 0
         for v in self.G.nodes:
+            # V1
+            to_add = set()
             for i in range(2,len(self.G.nodes)):
                 V = {w for w in self.G.neighbors(v)}
                 V.update([v])
@@ -146,8 +243,21 @@ class RootedConnectecKHopDominatingSet(ConnectedKHopDominatingSet):
                     if v in i_neighborhood:
                         min_ij_sep = ConnectedKHopDominatingSet.min_ij_separator(self.G, v, self.root, set(i_neighborhood))
                         if min_ij_sep:
-                            self.m.addConstr(gp.quicksum(self.nodes[s] for s in min_ij_sep) >= self.nodes[v])
+                            number_of_found_separators += 1
+                            # V1
+                            min_ij_sep_frozenset = frozenset(min_ij_sep)
+                            to_add.add(min_ij_sep_frozenset)
+                            
+                            # V2
+                            # self.m.addConstr(gp.quicksum(self.nodes[s] for s in min_ij_sep) >= self.nodes[v])
+            # V1
+            cummulated_to_add += len(to_add)
+            # V1
+            self.m.addConstrs(gp.quicksum(self.nodes[s] for s in min_ij_sep) >= self.nodes[v] for min_ij_sep in to_add)
+        
+        # print(cummulated_to_add, number_of_found_separators)
         
+    
     def add_simple_path_length_constraint(self):
         self.m.addConstrs((self.nodes[v] * len(nx.algorithms.shortest_path(self.G, self.root, v))) <= gp.quicksum(self.nodes) for v in self.G.nodes)
         # self.m.addConstrs((self.nodes[v] * self.nodes[w] * len(nx.algorithms.shortest_path(self.G, v, w))) <= gp.quicksum(self.nodes) for v in self.G.nodes for w in self.G.nodes)
@@ -162,10 +272,14 @@ class RootedConnectecKHopDominatingSet(ConnectedKHopDominatingSet):
 
     def intermediate_node_constraint(self, exclude):
         self.m.addConstrs(2 * self.nodes[v] <= gp.quicksum(self.nodes[w] for w in G.neighbors(v)) for v in G.nodes if v not in exclude)
+        
+    def getVars(self):
+        return self.m.getVars()
 
 
 if __name__ == '__main__':
     G = lp_to_nx_graph.read(sys.argv[1])
+    G = G.to_undirected()
     
     if(len(sys.argv) > 2):
         k = int(sys.argv[2])
@@ -173,7 +287,7 @@ if __name__ == '__main__':
         k = 1
         
     
-    # G = lp_to_nx_graph.read("/home/mario/Dokumente/Uni/Bachelorarbeit/git_repo/bachelor-mario-surlemont/python/lp graphs/asymmetric.lp")
-    # k = 2
-    dsProb = RootedConnectecKHopDominatingSet(G, k, 0)
+    # G = lp_to_nx_graph.read("/home/mario/Dokumente/Uni/Bachelorarbeit/git_repo/bachelor-mario-surlemont/python/lp graphs/small-leaf.lp")
+    # k = 1
+    dsProb = RootedConnectecKHopDominatingSet(G, k, 0, constraints = ConnectivityConstraint.SEPARATORS)
     dsProb.solve_and_draw()