diff --git a/kadmos/graph/graph_data.py b/kadmos/graph/graph_data.py
index a94828b4481b493e45f62ca3c77d98aadcb45f71..1fc0bff7f318b256e65842714e719754de94d187 100644
--- a/kadmos/graph/graph_data.py
+++ b/kadmos/graph/graph_data.py
@@ -358,6 +358,58 @@ class DataGraph(KadmosGraph):
 
         return
 
+    def get_coupling_matrix(self, function_order_method='manual', node_selection=None):
+        """Function to determine the role of the different functions in the FPG.
+
+        :param function_order_method: algorithm to be used for the order in which the functions are executed.
+        :type function_order_method: basestring
+        :param node_selection: selection of nodes for which the coupling matrix will be calculated only
+        :type node_selection: list
+        :return: graph with enriched function node attributes and function problem role dictionary
+        :rtype: FundamentalProblemGraph
+        """
+
+        # Make a copy of the graph, check it and remove all inputs and outputs
+        if node_selection:
+            graph = self.get_subgraph_by_function_nodes(node_selection)
+        else:
+            graph = self.cleancopy()
+        nodes_to_remove = list()
+        # TODO: Consider using the check function
+        assert not graph.find_all_nodes(subcategory='all problematic variables'), 'Graph still has problematic variables.'
+        nodes_to_remove.extend(graph.find_all_nodes(subcategory='all inputs'))
+        nodes_to_remove.extend(graph.find_all_nodes(subcategory='all outputs'))
+        graph.remove_nodes_from(nodes_to_remove)
+
+        # Determine and check function ordering method
+        assert function_order_method in self.OPTIONS_FUNCTION_ORDER_METHOD
+        if function_order_method == 'manual':
+            if node_selection:
+                function_order = node_selection
+            else:
+                assert 'function_order' in graph.graph['problem_formulation'], 'function_order must be given as attribute.'
+                function_order = graph.graph['problem_formulation']['function_order']
+        elif function_order_method == 'random':
+            function_order = graph.find_all_nodes(category='function')
+
+        # First store all the out- and in-edge variables per function
+        function_var_data = dict()
+        # noinspection PyUnboundLocalVariable
+        for function in function_order:
+            function_var_data[function] = dict(in_vars=set(), out_vars=set())
+            function_var_data[function]['in_vars'] = [edge[0] for edge in graph.in_edges(function)]
+            function_var_data[function]['out_vars'] = [edge[1] for edge in graph.out_edges(function)]
+        # Create an empty matrix
+        coupling_matrix = np.zeros((len(function_order), len(function_order)), dtype=np.int)
+        # Create the coupling matrix (including circular dependencies)
+        for idx1, function1 in enumerate(function_order):
+            for idx2, function2 in enumerate(function_order):
+                n_coupling_vars = len(set(function_var_data[function1]['out_vars']).
+                                      intersection(set(function_var_data[function2]['in_vars'])))
+                coupling_matrix[idx1, idx2] = n_coupling_vars
+
+        return coupling_matrix
+
     def get_possible_function_order(self, method, multi_start=None):
         """ Method to find a possible function order, in the order: pre-coupled, coupled, post-coupled functions
 
@@ -675,25 +727,16 @@ class DataGraph(KadmosGraph):
         # Calculate lower bound for each initial branch
         for node in nodes:
             node = [node]
-            lower_bound = self._get_lower_bound_branch_and_bound(node, nodes)
-            active_branches.append([node, lower_bound])
+            lower_bound_feedback, lower_bound_size = self._get_lower_bound_branch_and_bound(node, nodes)
+            active_branches.append([node, lower_bound_feedback, lower_bound_size])
 
         while True:
-            # Get lower bound information for each active branch
-            lower_bounds = np.array([branch[1] for branch in active_branches])
-
-            # Find indices of the branches with the lowest lower bound
-            min_lower_bound = np.min(lower_bounds)
-            indices_best_branches = [idx for idx, v in enumerate(lower_bounds) if v == min_lower_bound]
-
-            # Find the corresponding branches
-            best_branches = [active_branches[idx][0] for idx in indices_best_branches]
-
-            # Use a depth-first search to determine the branch that will be explored further
-            length_paths = [len(path) for path in best_branches]
-            idx_path = np.argmax(length_paths)
-            idx_selected_branch = indices_best_branches[idx_path]
-            selected_branch = active_branches[idx_selected_branch][0]
+            # Sort branches starting from low feedback and low size
+            sorted_branches = sorted(active_branches, key=lambda x: (x[1], x[2]))
+            # Find all branches with the least feedback and size
+            best_branches = [sorted_branches[idx][0] for idx in range(len(sorted_branches)) if sorted_branches[idx][1] == sorted_branches[0][1] and sorted_branches[idx][2] == sorted_branches[0][2]]
+            # Select the longest branch to further explore
+            selected_branch = max(best_branches, key=len)
 
             # Check whether the branch is a complete solution. If so, the best solution is found and iteration stopped
             if len(selected_branch) == len(nodes):
@@ -705,10 +748,12 @@ class DataGraph(KadmosGraph):
                 if node not in selected_branch:
                     node = [node]
                     current_order = selected_branch + node
-                    lower_bound = self._get_lower_bound_branch_and_bound(current_order, nodes)
-                    active_branches.append([current_order, lower_bound])
+                    lower_bound_feedback, lower_bound_size = self._get_lower_bound_branch_and_bound(current_order,
+                                                                                                    nodes)
+                    active_branches.append([current_order, lower_bound_feedback, lower_bound_size])
 
-            # Delete the explored branch
+            # Find index of explored branch in active branches and delete explored branch
+            idx_selected_branch = active_branches.index([selected_branch, sorted_branches[0][1], sorted_branches[0][2]])
             del active_branches[idx_selected_branch]
 
         function_order = selected_branch
@@ -734,39 +779,21 @@ class DataGraph(KadmosGraph):
             if node not in function_order:
                 function_order.append(node)
 
-        # Get coupling matrix for function order
-        # Make a subgraph of the nodes
-        subgraph = self.get_subgraph_by_function_nodes(function_order)
-        nodes_to_remove = []
-        nodes_to_remove.extend(subgraph.find_all_nodes(subcategory='all inputs'))
-        nodes_to_remove.extend(subgraph.find_all_nodes(subcategory='all outputs'))
-        subgraph.remove_nodes_from(nodes_to_remove)
-
-        # Get coupling matrix
-        # First store all the out- and in-edge variables per function
-        function_var_data = dict()
-        # noinspection PyUnboundLocalVariable
-        for function in function_order:
-            function_var_data[function] = dict(in_vars=set(), out_vars=set())
-            function_var_data[function]['in_vars'] = [edge[0] for edge in subgraph.in_edges(function)]
-            function_var_data[function]['out_vars'] = [edge[1] for edge in subgraph.out_edges(function)]
-        # Create an empty matrix
-        coupling_matrix = np.zeros((len(function_order), len(function_order)), dtype=np.int)
-        # Create the coupling matrix (including circular dependencies)
-        for idx1, function1 in enumerate(function_order):
-            for idx2, function2 in enumerate(function_order):
-                n_coupling_vars = len(set(function_var_data[function1]['out_vars']).
-                                      intersection(set(function_var_data[function2]['in_vars'])))
-                coupling_matrix[idx1, idx2] = n_coupling_vars
-
-        # Calculate lower bound
-        lowerbound = 0
-        for i in range(len(nodes)):
-            for j in range(len(branch)):
-                if i > j and coupling_matrix[i, j] != 0:
-                    lowerbound += 1
+        coupling_matrix = self.get_coupling_matrix(node_selection=function_order)
+
+        # Calculate lower bound for both feedback and size
+        feedback = 0
+        size = 0
+        for idx1 in range(len(nodes)):
+            for idx2 in range(len(branch)):
+                if idx1 > idx2 and coupling_matrix[idx1, idx2] != 0:
+                    feedback += coupling_matrix[idx1][idx2]
+                    if idx1 < len(branch):
+                        size += (idx1-idx2+1)*coupling_matrix[idx1][idx2]
+                    else:
+                        size += (len(branch)-idx2+1)*coupling_matrix[idx1][idx2]
 
-        return lowerbound
+        return feedback, size
 
     def get_feedback_info(self, function_order):
         """Function to determine the number of feedback loops for a given function order
@@ -777,38 +804,16 @@ class DataGraph(KadmosGraph):
         :rtype int
         """
 
-        # Make a subgraph of the nodes
-        subgraph = self.get_subgraph_by_function_nodes(function_order)
-        nodes_to_remove = []
-        nodes_to_remove.extend(subgraph.find_all_nodes(subcategory='all inputs'))
-        nodes_to_remove.extend(subgraph.find_all_nodes(subcategory='all outputs'))
-        subgraph.remove_nodes_from(nodes_to_remove)
-
-        # Get coupling matrix
-        # First store all the out- and in-edge variables per function
-        function_var_data = dict()
-        # noinspection PyUnboundLocalVariable
-        for function in function_order:
-            function_var_data[function] = dict(in_vars=set(), out_vars=set())
-            function_var_data[function]['in_vars'] = [edge[0] for edge in subgraph.in_edges(function)]
-            function_var_data[function]['out_vars'] = [edge[1] for edge in subgraph.out_edges(function)]
-        # Create an empty matrix
-        coupling_matrix = np.zeros((len(function_order), len(function_order)), dtype=np.int)
-        # Create the coupling matrix (including circular dependencies)
-        for idx1, function1 in enumerate(function_order):
-            for idx2, function2 in enumerate(function_order):
-                n_coupling_vars = len(set(function_var_data[function1]['out_vars']).
-                                      intersection(set(function_var_data[function2]['in_vars'])))
-                coupling_matrix[idx1, idx2] = n_coupling_vars
+        coupling_matrix = self.get_coupling_matrix(node_selection=function_order)
 
         # Determine number of feedback loops
         n_feedback_loops = 0
         n_disciplines_in_feedback = 0
-        for i in range(coupling_matrix.shape[0]):
-            for j in range(coupling_matrix.shape[0]):
-                if i > j and coupling_matrix[i, j] != 0:
-                    n_feedback_loops += 1
-                    n_disciplines_in_feedback += (i-j+1)
+        for idx1 in range(coupling_matrix.shape[0]):
+            for idx2 in range(coupling_matrix.shape[0]):
+                if idx1 > idx2 and coupling_matrix[idx1, idx2] != 0:
+                    n_feedback_loops += coupling_matrix[idx1][idx2]
+                    n_disciplines_in_feedback += (idx1-idx2+1)*coupling_matrix[idx1][idx2]
 
         return n_feedback_loops, n_disciplines_in_feedback
 
@@ -2221,58 +2226,6 @@ class FundamentalProblemGraph(DataGraph, KeChainMixin):
 
         return
 
-    def get_coupling_matrix(self, function_order_method='manual', node_selection=None):
-        """Function to determine the role of the different functions in the FPG.
-
-        :param function_order_method: algorithm to be used for the order in which the functions are executed.
-        :type function_order_method: basestring
-        :param node_selection: selection of nodes for which the coupling matrix will be calculated only
-        :type node_selection: list
-        :return: graph with enriched function node attributes and function problem role dictionary
-        :rtype: FundamentalProblemGraph
-        """
-
-        # Make a copy of the graph, check it and remove all inputs and outputs
-        if node_selection:
-            graph = self.get_subgraph_by_function_nodes(node_selection)
-        else:
-            graph = self.cleancopy()
-        nodes_to_remove = list()
-        # TODO: Consider using the check function
-        assert not graph.find_all_nodes(subcategory='all problematic variables'), 'Graph still has problematic variables.'
-        nodes_to_remove.extend(graph.find_all_nodes(subcategory='all inputs'))
-        nodes_to_remove.extend(graph.find_all_nodes(subcategory='all outputs'))
-        graph.remove_nodes_from(nodes_to_remove)
-
-        # Determine and check function ordering method
-        assert function_order_method in self.OPTIONS_FUNCTION_ORDER_METHOD
-        if function_order_method == 'manual':
-            if node_selection:
-                function_order = node_selection
-            else:
-                assert 'function_order' in graph.graph['problem_formulation'], 'function_order must be given as attribute.'
-                function_order = graph.graph['problem_formulation']['function_order']
-        elif function_order_method == 'random':
-            function_order = graph.find_all_nodes(category='function')
-
-        # First store all the out- and in-edge variables per function
-        function_var_data = dict()
-        # noinspection PyUnboundLocalVariable
-        for function in function_order:
-            function_var_data[function] = dict(in_vars=set(), out_vars=set())
-            function_var_data[function]['in_vars'] = [edge[0] for edge in graph.in_edges(function)]
-            function_var_data[function]['out_vars'] = [edge[1] for edge in graph.out_edges(function)]
-        # Create an empty matrix
-        coupling_matrix = np.zeros((len(function_order), len(function_order)), dtype=np.int)
-        # Create the coupling matrix (including circular dependencies)
-        for idx1, function1 in enumerate(function_order):
-            for idx2, function2 in enumerate(function_order):
-                n_coupling_vars = len(set(function_var_data[function1]['out_vars']).
-                                      intersection(set(function_var_data[function2]['in_vars'])))
-                coupling_matrix[idx1, idx2] = n_coupling_vars
-
-        return coupling_matrix
-
     def get_mg_function_ordering(self):
         """Method to determine the function ordering for MDAO graphs (FPG and MDG) based on an FPG.