Keywords: Branch-Price-and-Cut; Numerical Safety; Cycle Partitioning; Dynamic Programming; Bidirectional Search

1 Introduction↩︎

A cycle partition of an undirected graph \(G=(V,E)\) is a partition of the set of nodes \(V\) into disjoint cycles. In the length-constrained cycle partition problem (LCCP) introduced in [1], we are additionally given a critical time \(q_i > 0\) for every node \(i \in V\) and a travel time \(t_{i,j} \geq 0\) for every edge \(\{i,j\} \in E\). We call a cycle \(C=(i_0,i_1,\ldots,i_K=i_0)\) length-feasible if it satisfies the length constraint \[\label{eq:lc} t(C) := t_{i_0,i_1} + t_{i_1,i_2} + \ldots + t_{i_{K-1},i_K} \leq q(C) := \min \{ q_{i_0},q_{i_1},\ldots,q_{i_K} \},\tag{1}\] i.e., the total time to traverse the edges of the cycle must not exceed the critical time of any node in the cycle. In other words, an agent that continuously travels along the cycle will visit each node \(i\) of the cycle with frequency at least \(1/q_i\). The LCCP then amounts to computing the smallest number of disjoint, length-feasible cycles that cover all nodes.

The setting of LCCP is found to be useful, e.g., in applications where the objective is to find the smallest number of agents to conduct periodic surveillance or maintenance actions [2]. Similar types of length constraints can be found in the literature on minimizing cycle length in kidney exchange programmes [3], where the length is simply the number of nodes in the cycle. [1] show that the LCCP is NP-hard by reduction from the traveling salesman problem (TSP), and that no polynomial-time approximation algorithm exists. They present a compact mixed-integer programming (MIP) formulation based on the Miller, Tucker, and Zemlin (MTZ) [4] subtour-elimination constraints for the TSP. This compact model is able to solve instances with up to 29 nodes to proven optimality.

In a later work [5], the same authors introduced a MIP formulation based on exponentially many subtour elimination constraints (SEC), separated as cutting planes, together with valid inequalities derived from cliques in conflict hypergraphs, and an efficient primal heuristic. A branch-and-cut implementation of the SEC model based on a state-of-the-art commercial MIP solver was able to solve instances with up to 52 nodes. Both the MTZ and SEC models require additional auxiliary variables to count the number of used cycles in the objective function and to ensure that the total time of a cycle respects the critical time constraint of its nodes. Despite the mentioned improvements, the linear relaxations of these compact formulations remain weak. The formulations are moreover highly symmetric, which can hamper branching, although modern MIP solvers include dedicated symmetry-handling techniques that mitigate this to some extent.

The goal of our work is to overcome computational bottlenecks of the previous approaches by considering a set partitioning formulation based on cycle variables. Let \(\Omega\) denote the set of all length-feasible cycles, and let \(a_i^C\in \{0,1\}\) be a constant equal to \(1\) iff node \(i\) is contained in cycle \(C\). Introducing binary variables \(\lambda_C\) to indicate whether cycle \(C\in \Omega\) is used in the solution, we arrive at the set partitioning (SP) formulation \[\label{prob:imp} \begin{array}{llll} \min\;\; & \sum_{C\in \Omega} \lambda_C& \\\\ \text{s.t.} & \sum_{C\in \Omega} a_i^C\lambda_C&= 1 & \quad \text{for all }i \in V, \\\\ & \hfill\lambda_{C} &\in \{0,1\} & \quad \text{for all }C\in \Omega. \end{array} \tag{2}\] This reformulation of LCCP is inspired by formulations for the cardinality-constrained multi-cycle problem in kidney exchange [3], where a maximum-value packing of cycles is required, and by exact methods in many logistics applications such as vehicle routing, see, e.g., [6].

As in these applications, it is not viable to solve 2 explicitly due to the exponential number of variables. Instead, we apply column generation in order to solve the linear programming (LP) relaxation of 2 dynamically, and integrate this into an exact branch-and-price algorithm in order to solve LCCP instances to proven optimality. We refer to [7]–[10] for a general overview of this methodology.

A major computational task in this approach is the repeated solution of the so-called pricing problem for dynamically adding variables with negative reduced cost to the LP relaxation of 2 . Given travel times \(t_{i,j} \geq 0\) for \(\{i,j\}\in E\) and both critical times \(q_i > 0\) and weights \(\pi_i\in\mathbb{R}\) for \(i\in V\), the pricing problem amounts to computing a cycle \(C\) that satisfies the length constraint 1 and maximizes the sum of node weights. We call this the length-constrained prize-collecting cycle problem (LCPCCP). Note that problem data may be non-metric, i.e., travel times \(t_{i,j}\) may violate the triangle inequality, and that the node weights \(\pi_i\), which are derived as dual multipliers associated with the partitioning constraint of the node in the LP relaxation, are not restricted in sign. Even in the metric case and when critical times are assumed to be constant, LCPCCP is strongly NP-hard by reduction from the Hamiltonian cycle problem (1). The paper is organized as follows.

In 2 we present a new label-setting dynamic-programming algorithm for LCPCCP including dominance rules for reducing the search space and an efficient bidirectional search that allows us to enumerate cycles only half-way. In 3, we describe how to exploit the triangle inequality for locally metric input data, including per-node detection, pruning of zero-dual nodes, and preemptive feasibility checks. In 4, we adapt two families of cutting planes from the literature, subset row cuts and clique cuts, to the set partitioning formulation 2 of LCCP and discuss their management during column generation. In 5 we describe how we obtain numerically safe lower bounds without floating-point inaccuracies, and in 6 we integrate the resulting column generation scheme into a numerically safe branch-price-and-cut algorithm enhanced by symmetry breaking, heuristic pricing, parallel pricing, and early branching. In 7 we describe the results of our computational study to analyze the performance of the new approach on benchmark instances from the literature and to quantify the impact of the individual improvement techniques.

Our branch-price-and-cut algorithm is able to solve 52 instances with up to 76 nodes to proven optimality, closing \(14\) previously unsolved instances. Already its branch-and-price core is on average 14.7 times faster than the best previous approach (see 1), while in addition ensuring numerically safe bounds. Among the algorithmic contributions, the ingredients that exploit the specific structure of LCCP are, to the best of our knowledge, new: the symmetry of the undirected pricing problem, which lets bidirectional search extend labels in only one direction; symmetry breaking by critical time across the per-node pricing problems; and the exploitation of the triangle inequality for locally metric data.

A preliminary version of this work, introducing the set partitioning formulation and the branch-and-price algorithm, was presented at the INFORMS Optimization Society (IOS) Conference 2024 [11]. The present manuscript extends it with cutting-plane separation (subset-row and clique cuts) and a framework for numerically-safe dual bounds, together with a considerably expanded computational study.

2 Solving the LP Relaxation by Column Generation↩︎

Column generation is an advanced technique for solving large linear programs and is an essential subroutine in branch-and-price algorithms [7], [8]. The method has recently been found [12] to date back to [13], and was independently developed by [14].

In order to solve the LP relaxation of 2 by column generation, we first note that the upper bounds \(\lambda_C\leq 1\) are implied by the partitioning constraints and can be removed. This avoids introducing a separate dual variable for each bound, which would otherwise enter the reduced cost of every column. Without the upper bounds, the reduced cost depends only on the partitioning duals \(\pi_i\), which are unrestricted in sign. We call the resulting LP the master problem (MP). Second, we replace \(\Omega\) with a subset of columns \(\Omega' \subset \Omega\) to form the restricted master problem (RMP). Column generation then repeatedly optimizes (RMP), each time producing a primal-dual pair of solutions \((\lambda, \pi)\) and solving the so-called pricing problem to check whether cycles \(C\) with negative reduced cost \[\label{prob:pp} \textstyle\overline{c}(C) = 1 - \sum_{i \in V} a^C_i \pi_i = 1 - \sum_{i \in C} \pi_i\tag{3}\] exist in \(\Omega \setminus \Omega'\). If yes, a subset of these are added to \(\Omega'\) and (RMP) is reoptimized; if not, then \((\lambda, \pi)\) is optimal for (MP). We refer to [7], [8] for details on column generation. As explained above, the pricing problem of minimizing 3 amounts to LCPCCP, the length-constrained prize-collecting cycle problem.

2.1 Exact Pricing by Dynamic Programming↩︎

In the following, suppose we restrict LCPCCP by fixing one node \(s\in V\) to be contained in the cycle. Then the resulting problem closely resembles a resource-constrained shortest path problem (RCSPP) [16] typically solved in vehicle routing applications on a directed graph to find a path of minimum reduced cost that starts at and returns to a given node such that resources accumulated along arcs lie within a constant resource interval specified for each node. Similarly, we propose to solve LCPCCP for each fixed start node \(s\) using dynamic programming as an implicit enumeration scheme over all length-feasible cycles.

2.1.1 A Label-Setting Algorithm for LCPCCP↩︎

To this end, we form partial cycles (paths) and iteratively extend these paths with incident unvisited nodes until \(s\) is visited again, essentially closing the cycle. If at any point a partial cycle becomes infeasible (cannot be extended to a length-feasible cycle) or is proved to be dominated (would only lead to cycles with same or larger reduced cost), it is discarded, see 2.1.2. This process is repeated until all non-dominated length-feasible cycles are explored. If there are cycles with negative reduced cost, the ones with minimum reduced cost are returned.

Let \(P=(i_0,i_1, \dotsc, i_K)\) be a path in \(G=(V,E)\) with nodes \(i_0, \dotsc, i_K \in V\) and edges \(\{i_0,i_1\}, \dotsc, \{i_{K-1},i_K\} \in E\). The path \(P\) can be represented by a so-called label \(\ell= (\mathscr{N}, v, \overline{c}, t, q)\) where

• \(\mathscr{N}= \mathscr{N}(\ell) = \{i_1, \dotsc, i_K\}\) is the (unordered) set of visited nodes excluding the initial occurrence \(i_0=s\) of the start node (so a closed cycle, where \(i_K=s\), includes \(s\)),

• \(v = v(\ell) = i_K\) is the last node visited in \(P\),

• \(\overline{c}= \overline{c}(\ell) = 1 - \sum_{i \in \mathscr{N}} \pi_i\) is the reduced cost of \(P\),

• \(t = t(\ell) = \sum_{k=1}^K t_{i_{k-1},i_k}\) is the total travel time taken by \(P\), and

• \(q = q(\ell) = \min_{k=0,1, \dotsc, K} q_{i_k}\) is the minimum critical time of all nodes in \(P\).

The initial label for a starting node \(s\) is given as \(\ell_s=( \{\}, s, 1, 0, q_s)\). During the search, a label \(\ell= (\mathscr{N}, v, \overline{c}, t, q)\) can be extended using any edge \(\{v(\ell),j\}\) from the last node to a new node \(j \in V \setminus \mathscr{N}\). This creates a new label \(\ell^+\) where \(\mathscr{N}(\ell^+) = \mathscr{N}(\ell) \cup \{j\}\), \(v(\ell^+) = j\), \(\overline{c}(\ell^+) = \overline{c}(\ell) - \pi_j\), \(t(\ell^+) = t(\ell) + t_{v(\ell),j}\), and \(q(\ell^+) = \min\{q(\ell), q_j\}\). We call a label \(\ell^* \neq \ell_s\) fully-extended if it represents a cycle, i.e., if \(v(\ell^*) = s\). Note that LCCP allows singleton cycles (loops) as part of the solution; these are included here as label \((\{s\}, s, 1-\pi_s, 0, q_s)\). In order to recover the cycle corresponding to a fully-extended label, we also keep track of each label’s predecessor, i.e., the label from which it was extended.

2.1.2 Pruning the Search Space by Feasibility and Dominance↩︎

Besides the data structures used to represent and extend labels, the efficiency of the resulting dynamic programming method largely depends on how much of the search space needs to be explored in order to terminate with a provably optimal cycle. A label \(\ell\) can be discarded if it can be proven to only lead to length-infeasible cycles. This is the case if \(t(\ell) > q(\ell)\). Then for any extension \(\ell^+\) of label \(\ell\), we have \(t(\ell^+) \geq t(\ell) > q(\ell) \geq q(\ell^+)\). By induction, all fully-extended labels \(\ell^*\) obtained by successive extensions of label \(\ell\) satisfy \(t(\ell^*) > q(\ell^*)\), hence the corresponding cycles violate the length constraint 1 .

An effective technique for pruning the search space further is to exploit dominance relations. A label \(\ell_a\) is said to dominate label \(\ell_b\) if they have the same end node \(v=v(\ell_a)=v(\ell_b)\), and all cycles reachable from successive extensions of \(\ell_a\) have an equal or lower reduced cost than those reachable from successive extensions of \(\ell_b\). Then again, by induction, we may discard label \(\ell_b\) in the search, because at least one minimum reduced cost cycle will remain. Formally, this yields

The dominance rule can be strengthened by enriching the visited set with unreachable nodes. A node \(u\) is unreachable from a label \(\ell\) if visiting \(u\) would make it impossible to complete a feasible cycle, i.e., if \(t(\ell) + t_{v(\ell), u} > q(\ell)\). When extending a label to node \(j\), for each neighbor \(u\) of \(j\) that is unreachable from the extended label, we add \(u\) to the visited set. As these nodes cannot be part of any feasible extension, this does not affect correctness but strengthens dominance by making the visited set larger, allowing more labels to be pruned earlier in the search.

2.1.3 Bidirectional Search↩︎

Exploring the search space from both forward and backward directions simultaneously is another successful technique to reduce the number of explored labels. For the RCSPP, it has been shown that bidirectional search can be more efficient than monodirectional search [17], [18]. The efficiency gain comes from the fact that the number of non-dominated labels can grow exponentially with the length of the paths. Hence, if it is possible to explore forward and backward paths only until a halfway point and merge them to full paths, this can drastically reduce the total number of labels generated.

For LCPCCP, we can apply bidirectional search even more effectively than for the general RCSPP. Since the graph is undirected and the resource consumption by travel times is symmetric, only extensions from one direction are needed. Specifically, in bidirectional search for a fixed start node \(s\), we only extend labels \(\ell\) with \(t(\ell) < \frac{q(\ell)}{2}\). We then merge non-dominated labels as follows.

Let labels \(\ell_a\) and \(\ell_b\) both start at \(s\), end at \(v=v(\ell_a) = v(\ell_b)\), and be otherwise disjoint, i.e., \(\mathscr{N}(\ell_a) \cap \mathscr{N}(\ell_b) = \{v\}\). Then we create a cycle by concatenating the path represented by \(\ell_a\) and, in reverse order, the path represented by \(\ell_b\). We denote the newly merged label by \(\ell_a||\ell_b\), given by \(\mathscr{N}(\ell_a||\ell_b) = \mathscr{N}(\ell_a) \cup \mathscr{N}(\ell_b) \cup \{s\}\), \(v(\ell_a||\ell_b) = s\), \(\overline{c}(\ell_a||\ell_b) = \overline{c}(\ell_a) + \overline{c}(\ell_b) - \pi_s + \pi_{v(\ell_a)} - 1\), \(t(\ell_a||\ell_b) = t(\ell_a) + t(\ell_b)\), and \(q(\ell_a||\ell_b) = \min\{q(\ell_a), q(\ell_b)\}\), see 2 for an illustration. The reduced cost \(\overline{c}(\ell_a||\ell_b)\) is exactly the reduced cost \(1 - \sum_{i \in \mathscr{N}(\ell_a||\ell_b)} \pi_i\) of the merged cycle, where the three correction terms applied to the naive sum \(\overline{c}(\ell_a) + \overline{c}(\ell_b)\) account for the conventions in the label definition above. Each label’s reduced cost carries the constant \(1\) (the cost of one cycle), so the duplicate is removed by \(-1\). The shared end node \(v(\ell_a)=v(\ell_b)\) lies in both \(\mathscr{N}(\ell_a)\) and \(\mathscr{N}(\ell_b)\), so its dual is subtracted twice and added back once by \(+\pi_{v(\ell_a)}\). Finally, the start node \(s\) belongs to the closed merged cycle but, by convention, to neither \(\mathscr{N}(\ell_a)\) nor \(\mathscr{N}(\ell_b)\); its dual is therefore included by \(-\pi_s\). We only merge a disjoint pair whose resulting cycle is length-feasible, i.e., for which \(t(\ell_a) + t(\ell_b) \leq \min\{q(\ell_a), q(\ell_b)\}\).

We proceed to prove formally that considering the pairs of labels described above for merging is sufficient to obtain an optimal solution of LCPCCP, even in the presence of dominance checks.

Note that bidirectional search may produce more cycles than monodirectional search. The reason is that many of the labels at the halfway point could be dominated by labels found later in monodirectional search. Adding these potentially sub-optimal cycles to (RMP) can at the same time aid the convergence of column generation and make re-optimizing the LP much harder. We try to safeguard against this effect by sorting the cycles by reduced cost and returning at most a fixed number (\(50\) per starting node in our implementation) of the ones with lowest reduced cost.

2.2 Non-elementary Cycle Relaxation↩︎

The idea of solving a relaxation of the pricing problem in order to obtain valid lower bounds has been successfully applied in vehicle routing [19], [20]. The similarity of the pricing problem, RCSPP, to LCPCCP allows us to adapt this technique as follows. We call a cycle elementary if it visits every node at most once, and non-elementary if it revisits one or more nodes. We employ the so-called ng-paths relaxation [19], allowing some non-elementary cycles to be generated from the pricing problem. This gives us a lower bound on the minimum reduced cost, as it also contains the set of all elementary cycles.

In order to obtain a hierarchy of non-elementary relexations that converges to the result of LCPCCP with only elementary cycles, we define neighborhoods \(\mathcal{N}(i) \subseteq V\) for each node \(i \in V\), representing the “memory” of visited nodes maintained during label extension. When extending from node \(v\) to node \(j\), we record the visit (add \(j\) to the visited set) only if \(j \in \mathcal{N}(v)\) or \(v \in \mathcal{N}(j)\), and we forbid revisiting a node in the visited set. Hence a subcycle cannot occur between nodes that are in each other’s neighborhoods. Conversely, nodes outside mutual neighborhoods may be revisited, allowing non-elementary cycles. Setting \(\mathcal{N}(i) = V\) for all \(i\) recovers the elementary pricing algorithm. We initialize \(\mathcal{N}(i) = \{i, j^*\}\) where \(j^* = \arg\min_{j \neq i} t_{ij}\) is the closest neighbor of \(i\).

In the remainder of this section, we describe how subcycles in the generated cycles are detected, how the neighborhoods are expanded to prevent these subcycles from reoccurring, how a valid lower bound is obtained even while non-elementary cycles are generated, and how the relaxation interacts with bidirectional search.

Subcycle Detection.

When pricing returns a cycle \(C = (i_0, i_1, \ldots, i_K = i_0)\), we check whether it is elementary by detecting repeated nodes. We maintain a map from each node to its first occurrence index in the cycle. If node \(i_k\) was previously seen at index \(k' < k\), then \((i_{k'}, i_{k'+1}, \ldots, i_{k-1})\) forms a subcycle. All such subcycles are collected for neighborhood expansion.

Neighborhood Expansion.

For each detected subcycle \(S = (s_1, s_2, \ldots, s_m)\), we expand the neighborhoods of all involved nodes to include each other: \[\mathcal{N}(s_j) \leftarrow \mathcal{N}(s_j) \cup \{s_1, s_2, \ldots, s_m\} \quad \text{for all } j = 1, \ldots, m.\] This ensures that the same subcycle cannot reoccur: since all nodes in \(S\) are now mutually in each other’s neighborhoods, any future visit among them will be remembered and a revisit will be blocked.

After expanding neighborhoods, pricing is repeated. This process continues until all returned cycles are elementary, at which point the minimum reduced cost cycle is guaranteed to be elementary. The neighborhoods grow monotonically, and in the worst case converge to \(\mathcal{N}(i) = V\) for all \(i\), recovering the exact elementary algorithm.

Lagrangian Bound.

Importantly, even when non-elementary cycles are returned, the pricing problem provides a valid Lagrangian lower bound on the optimal LP objective. Let \(z^*_s\) denote the minimum reduced cost found for starting node \(s\) (possibly from a non-elementary cycle). Then \(z^*_{MP} \geq z^*_{RMP} + \sum_{s \in V} \min\{0, z^*_s\}\) is a valid lower bound on the master problem, since the relaxed pricing problem contains all elementary cycles.

Bidirectional Search.

The ng-relaxation is compatible with bidirectional search from 2.1.3 after relacing the merge condition in order to account for the partial tracking of the visited sets. In the elementary case, we require \(\tilde{\mathscr{N}}(\ell_a) \cap \tilde{\mathscr{N}}(\ell_b) = \{v\}\), where \(v = v(\ell_a) = v(\ell_b)\). However, with partially tracked visited sets, the meeting node \(v\) itself may not be tracked if its incident edges lie outside the mutual neighborhoods. We therefore relax the condition to \[\tilde{\mathscr{N}}(\ell_a) \cap \tilde{\mathscr{N}}(\ell_b) \subseteq \{v\},\] i.e., the two partial paths may share at most the meeting node \(v\), and possibly not even that if \(v\) is untracked. If they share any tracked node other than \(v\), a genuine overlap exists and the merge is rejected. Untracked nodes may still appear in both paths, potentially creating non-elementary cycles, which are handled by the neighborhood expansion procedure.

2.3 Greedy Pricing Heuristic↩︎

In order to avoid expensive labeling iterations, we first run a greedy heuristic to find improving cycles. The heuristic explores the neighborhood of each cycle \(C \in \Omega'\) currently in (RMP). For each cycle \(C\), we first remove all intermediate nodes with negative dual values, provided the resulting cycle remains length-feasible. Then, for each node \(i \notin C\) with positive dual \(\pi_i > 0\), we attempt to insert \(i\) at every position in the cycle, keeping all feasible insertions that yield negative reduced cost of the cycle.

This heuristic is run in parallel over all cycles in \(\Omega'\). The resulting improving cycles are sorted by reduced cost, and the best ones are added to (RMP). To avoid adding too many columns in a single iteration, we stop after adding at most \(|V|\) new columns. Only if no improving cycles are found do we proceed to the full labeling algorithm. We observed that this greedy approach is very effective at driving down the LP objective value early in the column generation process and significantly reduces the number of expensive labeling iterations required.

3 Exploiting the Triangle Inequality↩︎

The algorithm developed up to this point works for any choice of non-negative travel times. For many real-world applications, such as in aerial vehicle routing where distances are Euclidean, we know additionally that the triangle inequality holds, i.e., that \(t_{i,j} \leq t_{i,k} + t_{k,j}\) holds for all \(\{i,j\},\{i,k\},\{j,k\} \in E\). Under this assumption, the set partitioning formulation can be turned into a set covering formulation, which is easier to solve because the dual solution \(\pi\) becomes nonnegative. A solution for the set partitioning formulation can always be retrieved from a set covering solution by removing nodes that appear in more than one cycle from all but one cycle. This is guaranteed to be feasible since removing a node can only decrease the total travel time and increase the minimum critical time of the cycle.

Per-node detection.

Rather than assuming the triangle inequality holds globally, we detect for each node \(i \in V\) whether it can always be feasibly removed from any cycle. Specifically, node \(i\) satisfies the triangle inequality if \(t_{j,k} \leq t_{i,j} + t_{i,k}\) for all pairs \(j, k \in V\). Let \(V_{trieq} \subseteq V\) denote the set of such nodes. Note that a node enters \(V_{trieq}\) only when all required bypass edges exist (missing edges with travel time \(\infty\) fail the test), keeping it valid in general. For nodes in \(V_{trieq}\), we use set covering constraints (\(\geq 1\)) in the master problem 2 , while the remaining nodes retain set partitioning constraints (\(= 1\)). This mixed formulation allows us to exploit the triangle inequality even when it holds only partially. The effect of using covering constraints is that the dual variables of nodes in \(V_{trieq}\) are restricted to be nonnegative. This can weaken the LP bound but also aids convergence by stabilising the dual variables, preventing their oscillation in one direction.

Skipping zero-dual nodes.

By the arguments above, we can prove the following lemma showing that nodes with zero dual can be safely ignored during pricing.

Based on this lemma, we reduce the search space during label extension by skipping all neighbors \(j \in V_{trieq}\) with \(\pi_j = 0\).

Preemptive pruning.

Under the triangle inequality, we can also prune labels more aggressively. When extending a label to node \(j\), we check whether the cycle can still be closed feasibly by computing the time to directly return to the start node. If \(t(\ell) + t_{v(\ell), j} + t_{j, s} > q(\ell^+)\), we skip this extension, since no feasible cycle can be formed. This preemptive check avoids exploring labels that cannot lead to feasible cycles.

Post processing for set partitioning solutions.

The benefits above, namely the nonnegative and more stable duals and the pricing reductions of 2, all stem from solving with set covering constraints for the nodes in \(V_{trieq}\). Recovering a set partitioning solution from the covering solution afterwards is essentially free. Since a node in \(V_{trieq}\) may be covered by more than one cycle, we keep only its first occurrence and remove it from the remaining cycles. By the triangle-inequality argument given at the start of this section, each modified cycle stays feasible, and the number of cycles is unchanged; the result is therefore a valid set partitioning solution with the same objective, hence optimal for the original problem.

4 Cutting Planes↩︎

Given a fractional LP solution \(\lambda^*\), a cutting plane is a valid inequality that is violated by \(\lambda^*\) but satisfied by all integer feasible solutions. Strengthening the LP relaxation with the addition of valid inequalities has been shown to be effective [6], [21], [22]. In this section we present two popular families of cuts for set partitioning-based column generation, namely subset row cuts (SRCs) and clique cuts (CCs), followed by a discussion on whether generated cuts should be respected in the pricing routine.

4.1 Subset Row Cuts↩︎

Subset row cuts were introduced first by [21]. For a subset of nodes (rows) \(R \subseteq V\) (here representing cover constraint for each node) and \(k \in \{1, \dots, |R|\}\) , an SRC has the form \[\sum_{C\in \Omega} \lfloor \frac{1}{k} \sum_{i \in C} a_i^{C} \rfloor \lambda_{C} \leq \lfloor \frac{1}{k} |R| \rfloor\,.\] In our implementation we use \(k = 2\) and \(|R| = 3\). Intuitively, the cut is a Chvátal–Gomory rounding of the partitioning constraints of the rows in \(R\): summing the \(|R|\) constraints and dividing by \(k\) yields a fractional right-hand side that can be rounded down. For our choice \(k=2\) and \(|R|=3\), a cycle contributes to the left-hand side only if it covers at least two of the three nodes in \(R\). Since each of these nodes is covered exactly once in an integral partition, at most one selected cycle can cover two or more of them, which is precisely what the cut enforces; a fractional solution may instead spread the coverage of \(R\) across several cycles and thereby violate it. The separation routine iterates through all combinations of nodes and computes coefficients to find violated cuts.

To respect these generated cuts in the pricing problem, we augment each partial cycle \(\ell\) with one extra resource per active cut. Let \(S\) denote the set of currently active subset-row cuts, where each cut \(r \in S\) is generated by a node set \(R_r \subseteq V\). The added resource is a vector \[\rho(\ell) \in \{0, 1, \ldots, |R_r| - 1\}^{|S|}\] whose entry \(\rho_r(\ell) = |R_r \cap \mathscr{N}(\ell)|\) counts how many of the nodes in \(R_r\) the partial cycle has visited. We modify the accumulated reduced cost \(\overline{c}(\ell)\) to include the new duals by subtracting the corresponding dual \(\sigma_r\) for every \(k\) visits to the nodes in \(R_r\). The dominance rule is extended by the additional check \(\rho_r(\ell) \leq \rho_r(\ell')\) for all \(r \in S\), i.e., \(\ell\) must visit the nodes of each cut at most as many times as \(\ell'\) does in order to dominate it.

4.2 Clique Cuts↩︎

Clique cuts were first presented in [23] for the set partitioning problem and later used for a variant of the vehicle routing problem in [22], including how the dynamic programming pricing problem solver should be updated to respect the new duals from the generated cuts. Intuitively, two cycles that share a node can never be selected together in a partition, since that node would then be covered twice; hence among any set of cycles that pairwise conflict in this way, at most one may be used.

These cuts are based on cliques found on the so-called conflict graph. The conflict graph encodes which variables can cause infeasibility if both set to one in an integer solution. Each variable is represented by a node in this graph and an edge exists if setting both variables to 1 results in infeasibility. Consequently, every cycle variable is represented by a node, and an edge exists between a pair of nodes if the corresponding cycles visit any common node. A clique \(Q \subseteq \Omega\) on that graph represents a set of variables of which at most one can have value one in any feasible (integer) solution, giving the valid inequality \[\sum_{C\in Q} \lambda_C\leq 1\,.\]

We separate these cuts heuristically by enumerating \(4\)-combinations of the cycle variables in (RMP) and keeping those that form a clique in the conflict graph, i.e., every pair of cycles \(I, J\) in the combination shares at least one node, \(I \cap J \neq \emptyset\). We add all cuts that separate the current solution of (RMP).

4.3 Respecting Cuts in Pricing↩︎

Both SRCs and CCs are non-robust, i.e., they can remove solutions from the feasible region of the pricing problem. Robustness is not strictly required for guaranteeing optimality of 2 , since cuts (constraints) are not problem-defining, so the pricing problem can safely ignore them and still claim optimality [24]. However, respecting these cuts can lead to faster convergence of the column generation loop, with the tradeoff of making the pricing problem more complex. Therefore, it is not entirely obvious which approach is better.

Ignoring the generated cuts in the pricing problem allows us to generate many more cuts than would otherwise be possible before rendering the pricing problem intractable. This is the route we take for clique cuts: although they can be respected in the pricing problem as described in [22], doing so increases the complexity of the pricing problem, so we add them only to (RMP) and ignore their duals in pricing. For subset-row cuts, we additionally evaluate respecting them in the pricing problem (7.4).

Generally, even when a cut is ignored during pricing, it still constrains the newly generated columns in the master: as each column is added to (RMP), we lift every active cut to it by computing the column’s coefficient in that cut. This lifting recovers part of the cut’s strengthening effect at essentially no additional pricing cost.

4.4 Managing Cuts during Branch-and-Price↩︎

To balance the benefit of cutting planes against the overhead of separation and LP re-optimization, we limit the number of cuts added. In each separation round, we add at most \(10\) violated cuts per family (SRCs and CCs), and maintain a total limit of \(100\) cuts across the entire branch-and-bound tree. When the limit is reached, we rely on SCIP’s cut management to remove inactive cuts.

When branching creates child nodes, we inherit the cuts from the parent node. This ensures that the strengthened formulation is preserved throughout the tree without re-separating the same cuts. For newly generated columns added to (RMP), we compute and update their coefficients for all existing cuts, ensuring the cut constraints remain valid as the formulation evolves.

5 Numerically-safe Lowerbounds↩︎

Up until this point, we have assumed that (RMP) is solved exactly, i.e., without any numerical errors. However, in practice, (RMP) is typically solved with numerical tolerances using floating-point solvers, which can potentially lead to numerical errors through the accumulation of tiny roundoff errors. Ensuring numerical safety in LP-based branch-and-bound has been studied extensively [25]–[27], and also in the context of exact branch and price the need for numerically safe computation has been identified for capacitated vehicle routing [28] and bin packing [29]. In this section we follow the proposal made in the latter works to compute numerically safe bounds by scaling the dual values to integers and solving the pricing problem in exact integer arithmetic, and discuss how the DP-based pricing routine can be adapted to avoid errors arising from floating-point arithmetic.

Let us first consider the dual of (RMP), \[\label{prob:drmp} \begin{array}{llll} \max\;\; & \sum_{i \in V} \pi_i & \\[1ex] \text{s.t.} & \sum_{i \in C} \pi_i &\leq 1 & \quad \text{for all }C\in \Omega', \\[1ex] & \pi_i \in \mathbb{R} && \quad \text{for all }i \in V\,. \end{array}\tag{4}\] When solving (RMP) with a standard floating-point LP solver one would get a primal-dual pair of solutions \((\lambda, \pi)\) that satisfy primal and dual constraints within some tolerance. Concretely, we obtain an approximately feasible dual solution \(\pi\) of (DRMP) that satisfies \[\sum_{i \in C} \pi_i \leq 1 + \Delta \quad \text{for all }C\in \Omega'\] for some \(\Delta > 0\), typically between \(10^{-6}\) and \(10^{-9}\). In unsafe column generation, one would conclude optimality of the master problem when no columns with negative reduced cost exist, i.e., if all columns satisfy \(1 - \sum_{i \in C} \pi_i \geq -\Delta\), and use the dual objective \(\sum_{i \in V} \pi_i\) as a lower bound. However, if \(-\Delta \leq 1 - \sum_{i \in C} \pi_i < 0\), the branch-and-price result may be a suboptimal solution because (i) the lower bound is too large and (ii) improving columns may be missed.

We first focus on the first issue and discuss how to derive safe dual bounds given an infeasible dual solution \(\pi^{inf}\). The key idea is to repair this infeasibility by reducing the value of the left-hand side of the dual constraints. To ensure that no inaccuracies occur due to floating-point arithmetic, we first scale up and approximate the dual values by integers. We multiply the dual values by a large scaling factor \(M\) and round down, obtaining integers \[\label{equ:piint} \pi^{int} = \lfloor M \cdot \pi^{inf} \rfloor\,,\tag{5}\] then test dual feasibility exactly in integer arithmetic, i.e., whether \(\sum_{i \in C} \pi^{int}_i \leq M\) holds for all \(C\in \Omega'\). If the solution is infeasible, we reduce the scaling factor by a factor of ten, i.e., we set \(M \leftarrow M / 10\), and recompute 5 . This is repeated until an exactly dual feasible integer solution is found, or \(M < 1\), in which case the repair fails. Because the partitioning duals \(\pi_i\) are unrestricted in sign (see 4 ), the repair only needs to enforce the cycle constraints. When the triangle inequality is exploited via covering constraints for the nodes in \(V_{trieq}\) (3), the corresponding duals must additionally be nonnegative; we therefore change any negative covering dual to zero before scaling.

With this solution, the dynamic programming procedure in 2 can be run using integer arithmetic. If no negative reduced cost columns are found then one can conclude that the dual bound \(\sum_{i \in V} \pi^{int}_i/M\) obtained from the repaired dual solution is numerically safe. The procedure is summarized in 3.

Since this procedure progressively decreases the sum of the values of the dual solution, the safe dual bound is potentially worse than the unsafe dual bound. However, since the objective value is always integral, this would only lead to a different bound if the safe dual solution value crosses the integer boundary, i.e., if \(\lceil ( \sum_{i \in V} \pi_{i}^{int} ) / M \rceil < \lceil \sum_{i \in V} \pi_{i}^{inf} \rceil\).

It is worth mentioning that the cuts described in 4 are naturally numerically safe, since they are based solely on combinatorial arguments that are not affected by floating-point errors. Computing the cut violation of the current (RMP) solution \(\lambda\) in floating-point arithmetic can lead to incorrect determination whather a cut is violated. This may have the effect to not include violated cuts or include satisfied cuts and potentially affects performance, but does not compromise the correctness of the final result. Adding cuts leads to now dual variables in (DRMP). The exact dual feasibility check need not be modified to account for the new cut duals. Since the duals only appear with nonnegative coefficients in the dual constraints, and they have nonnegative values, a dual feasible solution can also be constructed from one that has positive values for the cut duals by setting these duals to zero.

Note that the difference between the safe and unsafe bound can have several effects. On the one hand, the safe bound may be weaker than the exact dual bound. If the LP relaxation at a node is integral and the unsafe bound already certifies optimality, a safe value that rounds up to one less than this integer may fail to prove optimality and force a branching. On the other hand, the unsafe bound may be invalidly strong. Since the dual solution is only feasible within tolerance, the unsafe value can slightly exceed the true optimum: a safe value of \(7.9999\) correctly rounds up to \(8\), whereas an unsafe value of \(8.0001\) would round up to \(9\), which is precisely the error that the exact repair prevents. In our experiments neither situation occurred: on every instance where both the safe and the unsafe configuration converge, the rounded safe and unsafe bounds agree (see 7.3), so ensuring numerical safety never changed a certified bound.

6 The Branch-and-Price Algorithm↩︎

Finally, let us discuss different aspects of how we extend the column generation method described in 2 into an exact branch-and-price algorithm.

Edge Branching.

By default, we use a standard branching strategy based on implicit edge variables as described in [6]. For a fractional LP solution, we define implicit edge variables \(x_{ij} = \sum_{C\in \Omega': \{i,j\} \in C} \lambda_C\) representing the total usage of edge \(\{i,j\}\) across all cycles. We branch on the most used fractional edge variable that does not belong to a singleton cycle. The rationale is that the most used edge would create the biggest disturbance in the subproblem, leading to an early fathoming due to infeasibility or finding a feasible integer solution.

In the down branch (\(x_{ij} = 0\)), we add edge \(\{i,j\}\) to a set of deleted edges and fix all cycle variables to zero that contain this edge. In the up branch (\(x_{ij} = 1\)), we add edge \(\{i,j\}\) to a set of fixed edges and fix all cycle variables to zero that contain node \(i\) or \(j\) but do not contain edge \(\{i,j\}\). The labeling algorithm must also respect these branching decisions. For deleted edges, we simply skip any extension along a deleted edge during label extension, and reject label merges in bidirectional search if the connecting edge is deleted. For fixed edges, we automatically extend labels along fixed edges: when a label reaches a node \(v\) with a fixed edge \(\{v, w\}\), we immediately extend to \(w\) without considering other neighbors. This ensures that all generated cycles respect the branching constraints while avoiding redundant exploration.

Ryan-Foster Branching.

Ryan-Foster [30] is an alternative branching technique for set partitioning master problems that branches on pairs of nodes rather than edges. For a pair of nodes \(i,j \in V\), we compute an implicit pair variable \(p_{ij} = \sum_{C\in \Omega_{ij}} \lambda_C^{*}\), where \(\Omega_{ij}\) is the set of cycles containing both \(i\) and \(j\). Given a fractional pair variable \(p_{ij}\), we branch to eliminate this fractionality.

We maintain sets \(M_i^n, C_i^n\) for each node \(i \in V\) at branch-and-bound node \(n\), representing nodes that must and cannot appear together with \(i\) in the same cycle, respectively. Let \(n'\) denote the parent of node \(n\):

• To enforce \(p_{ij} = 0\) (nodes in different cycles): fix all cycles containing both \(i\) and \(j\) to zero, and set \(C_i^n \leftarrow C_i^{n'} \cup \{j\}\), \(C_j^n \leftarrow C_j^{n'} \cup \{i\}\).

• To enforce \(p_{ij} = 1\) (nodes in same cycle): fix all cycles containing exactly one of \(i\) or \(j\) to zero, and set \(M_i^n \leftarrow M_i^{n'} \cup \{j\}\), \(M_j^n \leftarrow M_j^{n'} \cup \{i\}\).

The labeling algorithm must be modified to respect these branching constraints:

• We add resources \(m(\ell)\) and \(c(\ell)\) to each label, representing nodes that must and cannot be visited, respectively.

• When extending label \(\ell\) to node \(v\): \(m(\ell^+) = m(\ell) \cup M_v\) and \(c(\ell^+) = c(\ell) \cup C_v\).

• During extension we never visit a forbidden node, maintaining \(c(\ell) \cap \mathscr{N}(\ell) = \emptyset\); the requirement that all must-nodes are visited, \(m(\ell) \subseteq \mathscr{N}(\ell)\), is checked only when the cycle is closed.

• Dominance requires additional conditions: \(\ell'\) dominates \(\ell\) only if additionally \(m(\ell') \subseteq m(\ell)\) and \(c(\ell') \subseteq c(\ell)\).

• In bidirectional search, merging labels \(\ell_a\) and \(\ell_b\) additionally requires: \(c(\ell_a) \cap \mathscr{N}(\ell_b) = \emptyset\), \(c(\ell_b) \cap \mathscr{N}(\ell_a) = \emptyset\), and \(m(\ell_a) \cup m(\ell_b) \subseteq \mathscr{N}(\ell_a) \cup \mathscr{N}(\ell_b)\).

Node Selection and Processing.

For branch-and-bound node selection, we employ the best-estimate rule with plunging described in [31]. After selecting a branch-and-bound node, we solve the LP relaxation using column generation. In each pricing iteration, the dynamic programming algorithm from 2 is called from each starting node \(s\in V\). As these calls are independent, they are run in parallel.

Symmetry Breaking.

Since every cycle can be generated starting from any of its nodes, we apply symmetry breaking to avoid generating duplicate cycles from different pricing problems. For each call to the dynamic program from a fixed starting node \(s\), we restrict the search to only visit nodes with critical time at least as large as \(q_s\). For nodes with equal critical time, we use the node index as a tie-breaker, only allowing nodes \(j\) where either \(q_j > q_s\) holds or \(q_j = q_s\) and \(j > s\) hold.

This ensures that each cycle is generated exactly once: from the starting node with the smallest critical time (and smallest index among ties). Moreover, ordering by critical time rather than index balances the computational effort across pricing problems. The problems with smaller critical time \(q_s\) consider more nodes, but prune labels more agressively due to the length constraint. The problems with larger \(q_s\) prune labels less agressively, but in turn these pricing problems admit fewer nodes to start with, namely only those with critical time at least \(q_s\).

When branching decisions fix an edge \((s, j)\) to be in the solution, we can skip the entire pricing problem starting from \(s\) if \(j\) has smaller critical time (or equal critical time but smaller index), since any cycle using this edge would be found by the pricing problem starting from \(j\). Similarly, when using Ryan-Foster branching, if \(M_s\) contains a node with smaller critical time (or equal critical time but smaller index), the pricing problem starting from \(s\) can be skipped.

Pricing.

We start each pricing call with the greedy heuristic described in 2.3, returning any negative reduced cost columns found. If this fails to find a cycle with negative reduced cost, then a relaxed version of the pricing problem is run with limited neighborhoods that only include the closest node and the node itself (as described in 2.2). Whenever we get non-elementary cycles, we increase the neighborhood of each node to include nodes that appear together with that node in a subcycle. From this relaxation we always compute a Lagrangian lower bound [32]. Let \(z^*_s\) be the optimal objective value of the pricing problem for a fixed starting node \(s\), and let \(z^*_{MP}\), \(z^*_{RMP}\) be the optimal objective values of the master problem MP, i.e., the LP relaxation of 2 , and its restriction RMP, respectively, then \(z^*_{MP} \geq z^*_{RMP} + \sum_{s \in V} \min\{0, z^*_s\} =: LB_{lg}\). From integrality of the objective function, we can use \(\lceil LB_{lg} \rceil\) as a valid lower bound. Whenever no negative reduced cost cycle is found in this relaxation, we exit pricing and declare the LP as solved to optimality. The overall column generation procedure is summarized in 4. To improve convergence of the column generation loop, we employ dual stabilization using smoothed duals \(\tilde{\pi} = \alpha \pi + (1-\alpha) \pi^{prev}\), where \(\pi^{prev}\) is the dual solution from the previous iteration and \(\alpha \in (0,1]\) is a smoothing parameter.

Early Branching.

Another acceleration idea based on integrality, inspired by the technique in [33], is to skip pricing in some branch-and-bound node \(n\) by looking at the lower bound \(\lceil LB_{lg} \rceil\) of the parent node. Let \(z^*_{RMP}\) be the first RMP objective value for \(n\) obtained after removing columns due to branching. If \(\lceil z^*_{RMP} \rceil = \lceil LB_{lg} \rceil\), then computing \(z^*_{MP}\) exactly will not improve over the lower bound of the parent node. Therefore, we can skip pricing in this subproblem and perform early branching.

Farkas Pricing and RMP Initialization.

If the RMP at a subproblem becomes infeasible due to branching and removal of columns, we use Farkas pricing [34], [35] to generate new columns that render the RMP feasible again, or to prove that the MP is infeasible. The root node RMP is always feasible, because we initialize it with the trivially feasible singleton cycles and the cycles in the primal solution generated by the Most-Critical-Vertex-Based Heuristic (MCV) from [5].

Heuristic Branch-and-Price.

The exactness of the branch-and-price algorithm relies on solving each pricing problem to proven optimality. When relaxing this rule, we obtain a fast but inexact variant, which we use to compute strong primal solutions quickly. One possibility to solve the pricing problem heuristically is to relax the checks for dominance during labeling, thereby possibly pruning optimal LCPCPP solutions by dominance. Concretely, we drop the visited-set subset condition and the cut-resource conditions of 4 required for exact dominance and say that a label \(\ell'\) (heuristically) dominates \(\ell\) as soon as \(\overline{c}(\ell') \leq \overline{c}(\ell)\), \(t(\ell') \leq t(\ell)\) and \(q(\ell') \geq q(\ell)\). This discards far more labels and substantially accelerates pricing; in addition, the ng neighborhood-expansion step of 2.2 is skipped. Because pricing is no longer exact, the minimum reduced cost may be missed and the resulting Lagrangian bound is not valid. However, it can act as a strong primal heuristic to initialize an exact branch-and-price solve.

7 Computational Results↩︎

We implemented the branch-and-price algorithm described in 6 using the Rust wrapper of the branch-and-price framework SCIP [36]–[38]. The dynamic programming method from 2 is called from each starting node in parallel using the Rayon library [39]. The rest of SCIP is sequential. All experiments were run with a time limit of two hours on a cluster of identical machines, each equipped with two Intel(R) Xeon(R) Gold 6338 processors running at 2.0 GHz (up to 3.2 GHz), with 32 cores per socket for a total of 64 cores and 128 threads per node, and 1 TB of RAM. Each run had exclusive access to a node, and both our method and the bnc-sec baseline used up to 16 threads: our pricing through explicit parallelization and Gurobi through its default threading. We use the standard benchmark set from [5], which consists of 84 instances with 14 to 100 nodes.

Our experimental study is guided by the following questions:

1. How do the exact and heuristic branch-and-price variants compare against the branch-and-cut method of [5], both in terms of primal solution quality and dual bounds?

2. Which algorithmic components contribute the most to performance?

3. What is the cost and benefit of ensuring numerical safety?

4. Do cutting planes help to improve the branch-price-and-cut algorithm?

In order to answer these questions, in 7.1, we first evaluate overall performance, comparing the primal and dual sides of our approach against the branch-and-cut baseline. The best solutions found by the heuristic configurations in 7.1 serve as initial upper bounds for all subsequent experiments, in order to focus on the dual task of proving optimality. In 7.2, we analyze the impact of the algorithmic components, isolating the contribution of each through an ablation study and comparing alternative design choices. Our numerically safe implementation is compared against an unsafe floating-point version in 7.3. Finally, in 7.4 we investigate the effect of adding cutting planes.

Unless stated otherwise, all experiments use edge branching, see 6, and the numerically safe implementation described in 5. We report the shifted geometric mean (SGM) of solving times with a shift of 1 second, including time outs at the time limit. Within each table, the SGM is taken over the set of instances solved by at least one configuration listed in that table; consequently, absolute SGM values are not directly comparable across tables.

7.1 Overall Performance↩︎

We compare the primal and dual performance of the following configurations:

• bnc-sec: branch and cut with subtour elimination constraints from [5] based on Gurobi 13.0.1 [40]

• bnp-mcv: branch and price with bidirectional labeling (2.1.3), parallelization, symmetry breaking, and early branching (6), and greedy pricing (2.3), using only the MCV heuristic solutions from [5] as initial upper bounds

• heurbnp: branch and price with the heuristic pricing variant of 6, which relaxes the dominance rule to trade exactness for speed

• heurbnp-rf: same as heurbnp, but with Ryan-Foster branching

• bnp-full: same as bnp-mcv, but initialized with the best known solution across all heuristic configurations as initial upper bound

• bnc-full: same as bnc-sec, but warm-started with the same best known solution as bnp-full

The starting point for all methods is the greedy MCV heuristic from [5], which provides an initial feasible solution for each instance. The finally best known solution (BKS) for each instance is the best primal bound found by any of the heuristic configurations (MCV, heurbnp, heurbnp-rf). bnp-full is initialized with this BKS as the initial primal bound. Tab. 1 reports aggregated statistics on primal and dual performance of the five methods.

bnc-sec is able to improve the primal baseline given by the solution of the greedy MCV heuristic on 23 instances and bnp-mcv on 26 instances. In comparison, the heuristic branch-and-price method heurbnp proves particularly effective: by using a relaxed pricing algorithm that trades exactness for speed, it improves solutions on 41 instances with an average gap to the best known solution of only 2.1%. Adding Ryan-Foster branching (heurbnp-rf) improves solutions on 47 instances, six more than heurbnp, and lowers the average gap to the best known solution to 0.2%. Both heuristic variants are moreover fast: their shifted geometric mean running times over the instances solved by an exact configuration are 6.7 seconds (heurbnp) and 5.6 seconds (heurbnp-rf), below those of the exact methods, since they forgo proving optimality. Notably, neither heuristic dominates the other: heurbnp finds a strictly better solution than all other configurations on 1 instance, while heurbnp-rf does so on 8 instances, which motivates initializing bnp-full with the best solution found across all heuristic configurations.

On the dual side, bnp-mcv solves 46 instances to optimality compared to 37 for bnc-sec. This amounts to nine additional instances, and the branch-and-price variant bnp-mcv is moreover over 7 times faster than bnc-sec in terms of shifted geometric mean. Since both are initialized with the same MCV solution, this is a fair comparison. It is moreover conservative with respect to the solver: our branch-and-price runs on the academic solver SCIP, whereas bnc-sec uses the commercial solver Gurobi, which displays substantially faster performance on general mixed-integer programs.

The relaxed pricing in heurbnp also achieves the same (rounded) root dual bound as exact column generation on every instance where the latter converges. However, note that this heuristic “dual bound” comes with no formal guarantee of exactness. When the best known solutions from the heuristic configurations are fed back as initial upper bounds, bnp-full solves 51 instances to optimality, 14 more than bnc-sec, and is 14.7 times faster. To control for the stronger initial upper bound available to bnp-full, we also run the branch-and-cut baseline warm-started with the same best known solution, denoted bnc-full. Even with this advantage, bnc-full solves 40 instances, still 11 fewer than bnp-full, and remains over 10 times slower in shifted geometric mean. The largest instance solved by bnp-full features 76 nodes compared to 52 nodes for the branch-and-cut baseline.

To summarize, both primal and dual performance are significantly improved using branch and price. In the following experiments, we use bnp-full as the baseline configuration.

7.2 Analysis of Algorithmic Components↩︎

Our best version bnp-full combines the following techniques: bidirectional labeling (2.1.3), parallelization, symmetry breaking, and early branching (6), greedy pricing (2.3), and improved initial upper bounds from the heuristic branch-and-price variants (6). In order to quantify the performance impact of each component, we compare bnp-full against modified versions where one technique is disabled at a time:

• nosymbr: without symmetry breaking

• nobidir: with mono-directional labeling

• nopar: without parallel pricing

• noearly: without early branching

• nogreedy: without greedy pricing

• nobestsol: without improved initial upper bounds

• withtrieq: exploiting triangle inequality (3)

Tab. 2 reports the aggregated results including statistics on the root performance and reports the “virtual best” performance over all methods that could be achieved if we could choose the best method for each instance with perfect foresight. In the following, we first examine the quality of the root LP relaxation, which explains the overall behavior of the algorithm, then quantify the impact of the individual acceleration techniques, before taking a closer look at the pricing process and at three alternative design choices.

LP Relaxations.

The column generation loop of bnp-full converges to an optimal root LP solution for 52 instances, running into the memory or time limit in the remaining cases. The key advantage of the set partitioning formulation is evident in the tightness of its LP relaxation: 49 instances are solved at the root node without branching. Comparing the root dual bounds on the 52 instances where bnp-full converges, and measuring both against the same optimal objective, the set partitioning formulation attains an average gap of only 0.4%, compared to 30.3% for the SEC formulation of [5]. On 44 of these instances, bnp-full produces a strictly better root dual bound than bnc-sec. Furthermore, on all 51 instances solved to optimality, the root dual bound of the set partitioning formulation (rounded up due to integrality) is at least \(z^* - 1\), demonstrating the tightness of the LP relaxation. With most instances decided at the root, overall performance therefore hinges on how quickly the root LP relaxation is solved, i.e., on the efficiency of pricing, which is precisely what the following acceleration techniques target.

Impact of the Improvement Techniques.

Symmetry breaking proves to be the most impactful technique: disabling it results in a 3.6-fold slowdown and reduces the number of solved instances from 51 to 46. Bidirectional labeling is the second most effective technique, with monodirectional labeling resulting in a 3.2 times slowdown and 4 fewer instances solved. Parallelization provides a 2.2-fold speedup but only affects the number of solved instances marginally. Interestingly, early branching has negligible impact on these instances, with nearly identical solving times and the same number of instances solved. Disabling early branching (noearly) uniquely solves instance t84_brazil58, which bnp-full does not solve within the time limit. The virtual best of all ablation configurations solves 52 instances, one more than bnp-full.

Greedy Pricing Heuristic.

We next examine the two components that keep the individual pricing iterations cheap: the greedy pricing heuristic and the non-elementary relaxation. The greedy pricing heuristic described in 2.3 proves highly effective at finding negative reduced cost columns. Across all 84 instances, the greedy heuristic succeeds in 93% of pricing calls on average, contributing 77% of all columns added during column generation. The remaining 23% of columns are found by the labeling algorithm, which is invoked only when the greedy heuristic fails to find improving columns. Furthermore, the greedy heuristic reduces the number of neighborhood expansions in the ng-relaxation from an average of 7.5 per instance (without greedy) to 3.6 per instance (with greedy), indicating that greedy pricing often finds columns that would otherwise require expanding the ng-neighborhoods. Despite generating more total columns, bnp-full with greedy pricing is faster than without, as the greedy heuristic avoids the more expensive labeling algorithm in most pricing iterations. Disabling greedy pricing solves the same 51 instances with the same root bound quality but is about 13% slower, so greedy pricing improves convergence speed rather than the bounds themselves.

Non-elementary Relaxation.

The ng-relaxation described in 2.2 allows the pricing algorithm to generate non-elementary cycles with limited neighborhoods, expanding them only when necessary. Across all instances, neighborhoods grow to at most 67% of full size, with an average maximum of 38%. Notably, no instance requires neighborhoods to reach 100%, demonstrating that the ng-relaxation effectively reduces pricing complexity while still finding optimal solutions. On average, 7.5 neighborhood expansions are needed per instance when non-elementary cycles are detected. We do not include a configuration that disables the ng-relaxation in the ablation study, since disabling it amounts to forcing full (elementary) neighborhoods from the start, precisely the worst case the relaxation avoids. As the neighborhoods never grow beyond 67% of full size, this elementary variant would only ever be slower, so the relaxation is always beneficial.

Triangle Inequality.

Finally, we consider three alternative design choices (exploiting the triangle inequality, Ryan-Foster branching, and dual stabilization), each of which improves an isolated metric, yet fails to improve overall performance. Exploiting the triangle inequality (withtrieq) reduces root column generation iterations from 126 to 64 on average, yet only 49 instances are solved. This is because the set covering relaxation used for triangle inequality nodes can weaken the root dual bound: on 2 of the 52 instances where the root LP of both configurations converges, the dual bound is exactly 1 unit weaker. These are precisely the 2 instances that withtrieq fails to solve, and withtrieq does not solve any additional instances compared to bnp-full.

Branching Strategy.

In addition to the above experiments, we tested the effect of switching the branching rule to Ryan-Foster branching instead of edge branching (6). To expose enough branching to compare the two rules, we run both strategies without the improved initial upper bounds; with these bounds almost all instances are solved at the root, leaving too few branchings to observe. For this reason Ryan-Foster branching is reported here separately and is intentionally omitted from the ablation table above. Edge branching solves 46 instances (SGM 10.5s) while Ryan-Foster solves 45 (SGM 17.6s), making it 1.7 times slower. Over the 42 instances solved by both, 15 require branching: edge branching explores 8 nodes on average (max 61), while Ryan-Foster explores 14 nodes on average (max 99).

Dual Stabilization.

Furthermore, we tested dual stabilization to reduce the oscillation of dual values during column generation. While stabilization reduces the number of column generation iterations by 14% at the root node (from 126 to 108 on average) and by 18% overall, this barely affects solving times: the SGM is 6.3 seconds with stabilization versus 6.5 seconds without, and both configurations solve the same 51 instances.

7.3 Cost of Numerical Safety↩︎

Tab. 3 shows the impact of numerical safety on performance. The unsafe version bnp-unsafe is approximately 40% faster than the safe version bnp-full (4.5s vs 7.6s shifted geometric mean) and solves the same number of instances (51): the unsafe version closes t62_hk48, while the safe version instead closes at62_p43. No dual bound failures were detected: on all instances where both versions converge, the safe and unsafe dual bounds agree. This demonstrates that numerical safety can be ensured with a moderate overhead, preserving the substantial performance gains of the branch-and-price approach, but that it may be beneficial to run the heuristic versions heurbnp and heurbnp-rf in unsafe mode.

7.4 Branch-Price-and-Cut Performance↩︎

In order to evaluate the options described in 4, we compare five different branch-price-and-cut configurations:

• bnp-full: no cutting planes (baseline)

• withcc: clique cuts without change to the pricing problem

• withsrc: subset row cuts without change to the pricing problem

• withsrc+cc: both subset row and clique cuts without change to the pricing problem

• withsrc+pricing: subset row cuts with modifications to the pricing problem

Tab. 4 reports aggregated results including statistics on the root performance and again reports the “virtual best” performance over these methods.

First, cutting planes are separated on very few instances: subset-row cuts are added on only 3 of the 84 instances (t62_dantzig42, t84_eil51 and t84_eil76, with 11 cuts in a single separation round each), and clique cuts on 2 instances (t84_eil51 and t84_eil76, with 33 and 26 cuts over 5 and 4 rounds, respectively). This is a direct consequence of the tight LP relaxation: most instances are solved at the root with little fractionality left to separate. Adding cutting planes to the RMP without modifying the pricing problem generally fails to help: clique cuts alone (withcc) are roughly neutral, even marginally faster (a 3% speedup), while solving the same instances, whereas subset row cuts (withsrc) cause a 27% slowdown and combining both (withsrc+cc) also a 27% slowdown. Adding cuts only to the RMP can even be counterproductive: relative to bnp-full, withsrc fails to solve t62_att48, t62_dantzig42 and t84_eil51 (48 solved), while withsrc+cc fails on t62_dantzig42 and t84_eil51 (49 solved). The number of column generation iterations remains nearly identical across configurations on instances solved by all methods, so the slowdown is attributable to the overhead of cut separation and LP re-optimization, which is not compensated by tighter bounds.

Respecting the subset row cuts in the pricing problem (withsrc+pricing) reverses this effect on t84_eil51: the cuts close the root node, so the instance is solved without branching (in 13.7 s, versus the 25.5 s and 15 branch-and-bound nodes needed by bnp-full), recovering the instance that withsrc loses. By generating columns that already satisfy the cuts, this configuration solves 52 instances, one more than bnp-full (it additionally closes t62_hk48), with a 4% lower shifted geometric mean time (7.2 versus 7.6 seconds), the best of all cutting-plane configurations. The virtual best over all cutting plane configurations solves 52 instances, the same as withsrc+pricing.

8 Conclusion↩︎

To summarize, our computational results show that a branch-price-and-cut approach significantly improves over the previous state-of-the-art in solving the length-constrained cycle partition problem, both in terms of primal solution quality and dual performance. The best-performing branch-price-and-cut algorithm (withsrc+pricing) solves 52 instances from the standard benchmark set, compared to 40 for the previous branch-and-cut method, and scales to instances with up to 76 nodes (vs. 52 previously reported). The main driver of this performance improvement is the set partitioning formulation: it provides substantially tighter LP relaxations, with an average root gap of 0.4% compared to 30.3% for the subtour elimination formulation.

As a result, most instances are solved at the root node without branching. A natural question is whether the set-partitioning formulation satisfies the round-up property, i.e., whether the optimal integer objective always equals the LP relaxation value rounded up, \(z^* = \lceil z_{LP} \rceil\). Our results show that this does not hold in general: on t84_eil51 the LP value rounds up to 8 while the integer optimum is 9, which is exactly why branching is required there. On the other 50 of the 52 instances whose root LP converges the property does hold, and on the single open instance t84_eil76 it cannot be decided, as its LP rounds up to 8 while the best known solution uses 9 cycles. This near-universal round-up behavior is nonetheless what allows most instances to be solved at the root once a matching primal solution is found. For a study of the conditions under which the round-up property holds, see [41].

Among the algorithmic components, symmetry breaking, bidirectional labeling, and improved initial upper bounds contribute the most. The greedy pricing heuristic succeeds in 93% of pricing calls and produces root LP bounds matching the exact algorithm, while avoiding the more expensive labeling procedure. Subset row cuts respected in pricing are the only cutting plane configuration that improves performance, reducing the solving time without losing any instances. When branching is needed, edge branching proves more robust than Ryan-Foster branching, producing smaller search trees. Interestingly, techniques that reduce the number of column generation iterations, namely dual stabilization and triangle inequality exploitation, do not improve solving times, suggesting that the pricing problem is solved efficiently enough that fewer iterations do not compensate for the added overhead.

All results are based on a numerically safe combination of floating-point and integer arithmetic. While we did not notice any incorrect results when using unsafe floating-point computation, we could confirm that numerical safe computation incurs only a moderate overhead and preserves the substantial performance gains of the branch-price-and-cut approach.

The most pressing issue for future work is the efficiency of the pricing routine: instances beyond 76 nodes remain unsolved, and pricing dominates the running time. Closing these instances will likely require stronger dual bounds or a more scalable pricing algorithm. Recent advances in accelerating dynamic-programming-based pricing, such as parallel labeling for the resource-constrained shortest path problem [42] or GPU-accelerated state expansion [43], offer a promising avenue. The potential of further improvements, such as parallelizing the remaining sequential parts of the branch-and-bound search or separating additional cutting planes, can only be unlocked once pricing is no longer the bottleneck.

Funding↩︎

Research reported in this paper was partially supported through the Research Campus MODAL funded by the German Federal Ministry of Education and Research (fund numbers 05M14ZAM, 05M20ZBM) and the Deutsche Forschungsgemeinschaft (DFG) through the DFG Cluster of Excellence MATH+.

Data availability↩︎

The benchmark instances and the scripts used to generate all reported results are openly available in the lccp-ejor repository at https://github.com/mmghannam/lccp-ejor.

9 Per-Instance Results↩︎

Table 5: Per-instance results comparing bnp-full (branch-price-and-cut) with bnc-sec (branch-and-cut with subtour elimination). BKS: best known solution; PB: primal bound; DB: dual bound; DBr: dual bound at root; Gap: optimality gap (%); Time: solving time in seconds; Nodes: branch-and-bound nodes; Status: opt (optimal) or TL (time limit).

			bnp-full							bnc-sec
Instance	\(n\)	BKS	PB	DBr	DB	Gap	Time	Nodes	St.	PB	DB	Gap	Time	Nodes	St.
Table – continued from previous page
			bnp-full							bnc-sec
Instance	\(n\)	BKS	PB	DBr	DB	Gap	Time	Nodes	St.	PB	DB	Gap	Time	Nodes	St.
Continued on next page
t62_burma14	14	5	5	5	5	0	0.1	1	opt	5	5	0	0.1	43	opt
t84_burma14	14	6	6	6	6	0	0.1	1	opt	6	6	0	0.0	1	opt
t62_ulysses16	16	4	4	4	4	0	0.1	1	opt	4	4	0	0.0	1	opt
t84_ulysses16	16	6	6	6	6	0	0.1	1	opt	6	6	0	0.1	1	opt
at62_br17	17	5	5	5	5	0	0.1	1	opt	5	5	0	0.1	1	opt
at84_br17	17	6	6	6	6	0	0.0	1	opt	6	6	0	0.0	1	opt
t62_gr17	17	5	5	5	5	0	0.1	1	opt	5	5	0	0.2	1	opt
t84_gr17	17	8	8	8	8	0	0.1	1	opt	8	8	0	0.0	1	opt
t62_gr21	21	5	5	5	5	0	0.1	1	opt	5	5	0	0.8	34	opt
t84_gr21	21	8	8	8	8	0	0.1	1	opt	8	8	0	0.5	21	opt
t62_ulysses22	22	5	5	5	5	0	0.3	1	opt	5	5	0	0.1	1	opt
t84_ulysses22	22	7	7	7	7	0	0.1	1	opt	7	7	0	0.7	67	opt
t62_gr24	24	5	5	5	5	0	0.5	1	opt	5	5	0	1.1	25	opt
t84_gr24	24	7	7	7	7	0	0.1	1	opt	7	7	0	0.2	1	opt
t62_fri26	26	6	6	6	6	0	1.0	1	opt	6	6	0	1.9	158	opt
t84_fri26	26	8	8	8	8	0	0.1	1	opt	8	8	0	1.1	154	opt
t62_bayg29	29	5	5	5	5	0	1.1	1	opt	5	5	0	12.7	8494	opt
t62_bays29	29	6	6	6	6	0	0.6	1	opt	6	6	0	3.4	814	opt
t84_bayg29	29	8	8	8	8	0	0.1	1	opt	8	8	0	17.1	5969	opt
t84_bays29	29	8	8	8	8	0	0.1	1	opt	8	8	0	7.4	1018	opt
at62_ftv33	33	7	7	7	7	0	1.2	1	opt	7	7	0	3113.6	3353943	opt
at84_ftv33	33	9	9	9	9	0	0.1	1	opt	9	9	0	13.0	403	opt
at62_ftv35	35	6	6	6	6	0	3.8	1	opt	6	6	0	3285.3	14828814	opt
at84_ftv35	35	9	9	9	9	0	0.3	1	opt	9	9	0	165.6	26250	opt
at62_ftv38	38	6	6	6	6	0	9.6	1	opt	6	5	16.7	TL	36061643	TL
at84_ftv38	38	9	9	9	9	0	0.7	1	opt	9	9	0	546.9	278079	opt
t62_dantzig42	42	6	6	6	6	0	113.4	41	opt	6	6	0	1024.9	2590519	opt
t62_swiss42	42	6	6	6	6	0	87.5	1	opt	6	6	0	1630.6	3776539	opt
t84_dantzig42	42	9	9	9	9	0	0.9	1	opt	9	9	0	1364.0	253146	opt
t84_swiss42	42	9	9	9	9	0	0.9	1	opt	9	9	0	233.1	53847	opt
at62_p43	43	1	1	1	1	0	297.3	1	opt	1	1	0	0.1	0	opt
at84_p43	43	2	2	2	2	0	3.0	1	opt	2	2	0	0.1	1	opt
at62_ftv44	44	6	6	6	6	0	77.4	1	opt	7	5	28.6	TL	10118677	TL
at84_ftv44	44	9	9	9	9	0	0.7	1	opt	9	9	0	6617.2	1007542	opt
at62_ftv47	47	6	6	6	6	0	215.6	1	opt	7	6	14.3	TL	6819795	TL
at84_ftv47	47	8	8	8	8	0	1.2	1	opt	8	8	0	5322.9	6130397	opt
at62_ry48p	48	6	-	-	-	100.0	TL	-	TL	6	5	16.7	TL	11727846	TL
at84_ry48p	48	8	8	8	8	0	7.6	1	opt	9	8	11.1	TL	6659592	TL
t62_att48	48	6	6	6	6	0	4599.6	1	opt	6	6	0	1096.1	2278842	opt
t62_gr48	48	6	6	6	6	0	98.2	1	opt	7	5	28.6	TL	14151852	TL
t62_hk48	48	6	-	-	-	100.0	TL	-	TL	6	6	0	6026.9	20699967	opt
t84_att48	48	8	8	8	8	0	26.7	1	opt	8	8	0	1010.4	186278	opt
t84_gr48	48	8	8	8	8	0	5.6	1	opt	9	8	11.1	TL	9427085	TL
t84_hk48	48	9	9	9	9	0	3.6	1	opt	9	9	0	1904.9	540077	opt
t62_eil51	51	6	6	6	6	0	947.2	1	opt	7	5	28.6	TL	13281898	TL
t84_eil51	51	9	9	8	9	0	25.5	15	opt	9	8	11.1	TL	9935924	TL
t62_berlin52	52	6	-	-	-	100.0	TL	-	TL	7	6	14.3	TL	14358987	TL
t84_berlin52	52	9	-	-	-	100.0	TL	-	TL	9	8	11.1	TL	5048215	TL
at62_ft53	53	5	-	-	-	100.0	TL	-	TL	7	4	42.9	TL	8070726	TL
at84_ft53	53	8	8	8	8	0	47.4	1	opt	10	6	40.0	TL	2166224	TL
at62_ftv55	55	6	6	6	6	0	448.8	1	opt	7	5	28.6	TL	7753333	TL
at84_ftv55	55	9	9	9	9	0	9.2	1	opt	11	7	36.4	TL	648741	TL
t62_brazil58	58	5	-	-	-	100.0	TL	-	TL	6	5	16.7	TL	15534631	TL
t84_brazil58	58	8	-	-	-	100.0	TL	-	TL	8	7	12.5	TL	8560266	TL
at62_ftv64	64	7	-	-	-	100.0	TL	-	TL	8	5	37.5	TL	2958768	TL
at84_ftv64	64	10	10	10	10	0	26.2	1	opt	11	7	36.4	TL	358873	TL
at62_ft70	70	6	-	-	-	100.0	TL	-	TL	6	4	33.3	TL	4094939	TL
at62_ftv70	70	9	-	-	-	100.0	TL	-	TL	9	5	44.4	TL	2445597	TL
at84_ft70	70	9	-	-	-	100.0	TL	-	TL	9	6	33.3	TL	1727481	TL
at84_ftv70	70	9	9	9	9	0	177.6	1	opt	11	7	36.4	TL	872838	TL
t62_st70	70	6	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	6432455	TL
t84_st70	70	9	9	9	9	0	76.2	1	opt	9	7	22.2	TL	401218	TL
t62_eil76	76	6	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	3945851	TL
t62_pr76	76	6	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	3701310	TL
t84_eil76	76	9	9	8	8	11.1	TL	5	TL	10	7	30.0	TL	730776	TL
t84_pr76	76	9	9	9	9	0	1412.8	1	opt	10	7	30.0	TL	137363	TL
t62_gr96	96	7	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	1917562	TL
t84_gr96	96	9	-	-	-	100.0	TL	-	TL	10	6	40.0	TL	59244	TL
t62_rat99	99	7	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	1430164	TL
t84_rat99	99	10	-	-	-	100.0	TL	-	TL	10	7	30.0	TL	146129	TL
t62_kroA100	100	8	-	-	-	100.0	TL	-	TL	8	5	37.5	TL	426410	TL
t62_kroB100	100	8	-	-	-	100.0	TL	-	TL	8	5	37.5	TL	618315	TL
t62_kroC100	100	7	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	461929	TL
t62_kroD100	100	7	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	1541628	TL
t62_kroE100	100	7	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	1361385	TL
t62_rd100	100	7	-	-	-	100.0	TL	-	TL	7	5	28.6	TL	369576	TL
t84_kroA100	100	9	-	-	-	100.0	TL	-	TL	10	6	40.0	TL	16871	TL
t84_kroB100	100	9	-	-	-	100.0	TL	-	TL	10	6	40.0	TL	4821	TL
t84_kroC100	100	9	-	-	-	100.0	TL	-	TL	11	6	45.5	TL	6357	TL
t84_kroD100	100	9	-	-	-	100.0	TL	-	TL	10	7	30.0	TL	24214	TL
t84_kroE100	100	9	-	-	-	100.0	TL	-	TL	11	6	45.5	TL	17521	TL
t84_rd100	100	10	-	-	-	100.0	TL	-	TL	11	6	45.5	TL	4416	TL
at62_kro124p	124	7	-	-	-	100.0	TL	-	TL	7	4	42.9	TL	202101	TL
at84_kro124p	124	9	-	-	-	100.0	TL	-	TL	10	6	40.0	TL	71479	TL

References↩︎

---

1. ghannam@zib.de↩︎

2. gleixner@htw-berlin.de↩︎

3. edward.lam@monash.edu↩︎

4. mexi@zib.de↩︎