Integral column generation for the set partitioning problem

The integral simplex using decomposition (ISUD) algorithm was recently developed to solve efficiently set partitioning problems containing a number of variables that can all be enumerated a priori. This primal algorithm generates a sequence of integer solutions with decreasing costs, leading to an optimal or near-optimal solution depending on the stopping criterion used. In this paper, we develop an integral column generation (ICG) heuristic that combines ISUD and column generation to solve set partitioning problems with a very large number of variables. Computational experiments on instances of the public transit vehicle and crew scheduling problem and of the airline crew pairing problem involving up to 2000 constraints show that ICG clearly outperforms two popular column generation heuristics (the restricted master heuristic and the diving heuristic). ICG can yield optimal or near-optimal solutions in less than 1 hour of computational time, generating up to 300 integer solutions during the solution process.

Authors and Affiliations

1. GERAD and Department of Mathematics and Industrial Engineering, Polytechnique Montréal, Montréal, Canada

   Adil Tahir, Guy Desaulniers & Issmail El Hallaoui

Corresponding author

Correspondence to Adil Tahir.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix A: compatibility matrices

Variable substitution can be used to reduce the density of the constraint coefficient matrix of the CP (7)–(11). A first reformulation, proposed by Zaghrouti et al. (2014), is obtained by isolating each \(\lambda _l\) variable in the last constraint (8) in which it appears before substituting it in the rest of the formulation. This substitution yields the following reformulation of the CP:

where, for \(U\subset N\), \({\mathbf {c}}_U\) is the subvector of the cost coefficient vector in (7) associated with the variables \(v_j\), \(j\in U\). Constraint (17) can be written in the form \(({\mathbf {M}}_1{\mathbf {A}}_{I_S})^{\top }{\mathbf {v}}= 0\), where \({\mathbf {M}}_1\) is a compatibility matrix according to the following definition.

Definition 3

A matrix \({\mathbf {M}}\) is said to be a compatibility matrix if and only if \({\mathbf {M}}{\mathbf {A}}_j = 0\) for all compatible columns \({\mathbf {A}}_j\), \(j\in C_S\), and \({\mathbf {M}}{\mathbf {A}}_j \ne 0\) for all incompatible columns \({\mathbf {A}}_j\), \(j\in I_S\).

There exists an infinite number of compatibility matrices that can be used to define the CP. The choice of this matrix may have a significant impact on the computational times. We have chosen to use the matrix \({\mathbf {M}}_2\) introduced by Bouarab et al. (2017) that is specialized for vehicle routing and crew-scheduling problems. In these problems, each column is associated with a route or a crew schedule that performs a subset of tasks in a given sequence. The \({\mathbf {M}}_2\) matrix allows to measure the deviation of the incompatible columns with respect to the task sequences defined by the columns in solution S. This matrix is such that, if an incompatible column \({\mathbf {A}}_j\), \(j\in I_S\), does not respect the task ordering in a sequence, \({\mathbf {M}}_2{\mathbf {A}}_j\) contains a component equal to -1 each time that \({\mathbf {A}}_j\) covers a task but not its predecessor and a component equal to 1 each time that \({\mathbf {A}}_j\) covers a task but not its successor (see Bouarab et al. 2017 for details).

Appendix B: multi-phase strategy

Inspired from that developed by Elhallaoui et al. (2010) in the context of DCA, Zaghrouti et al. (2014) have elaborated a multi-phase acceleration strategy for ISUD which also aims at reducing the density of the CP constraint coefficient matrix. Each time that the CP needs to be solved, a sequence of phases can be invoked. Each phase is defined by a parameter k (a positive integer) which restricts the CP to a subset \(I^k_S\subseteq I_S\) of the incompatible columns such that \(I^k_S \subseteq I^{\ell }_S\) if \(k < \ell\). The sequence of phases is predetermined and corresponds to an increasing sequence \(k_1, k_2, \ldots , k_p\) of values of k with \(k_p = \infty\) (for instance, \(k = 2, 3, 4, 5, 6, \infty\)), where phase \(k_p = \infty\) means that \(I^{\infty }_S = I_S\). Starting in phase \(k_1\), the CP restricted to the subset \(I^{k_1}_S\) is solved. If its optimal value is negative, the process stops and the computed linear combination of incompatible columns is returned. Otherwise, the next phase is invoked. The process repeats until obtaining a negative optimal value for the CP or reaching phase \(k_p = \infty\).

The definition of a subset \(I^k_S\) used in phase k is based on a distance between an incompatible column and the vector subspace generated by the columns in the solution. This distance, called the degree of incompatibility, is defined as follows.

Definition 4

The degree of incompatibility\(\delta _j\) of an incompatible column \({\mathbf {A}}_j\), \(j\in I_S\), is given by \(\delta _j = || {\mathbf {M}}{\mathbf {A}}_j ||\), where \({\mathbf {M}}\) is a compatibility matrix (as defined in Appendix A).

In phase k of the multi-phase strategy for solving the CP, the index subset of the incompatible columns considered in the CP is defined by \(I^k_S = \{ j \in I_S \, | \, \delta _j \le k \}\). Bouarab et al. (2017) has proven the following result.

Proposition 2

In phase k , each incompatible column \({\mathbf {A}}_j\) , \(j\in I^k_S\) , has at most \(k+1\) non-zero coefficients in the constraint coefficient matrix of the CP if \({\mathbf {M}}_2\) is used as the compatibility matrix.

Consequently, for phases with a low k value, the constraint coefficient matrix has a low density. This helps to reduce the computational effort for solving the CP and also to produce integer directions without branching. Finally, it may allow to consider a very large number of incompatible columns in the CP if many are available.

Appendix C: computation of a complete dual solution

To compute a complete dual solution \({\varvec{\alpha }}= \left( \alpha _i \right) _{i\in T} \in {\mathbb {R}}^m\) to the MP, we use as in Bouarab et al. (2017) the dual variable values \({\hat{{\varvec{\pi }}}} \in {\mathbb {R}}^{m-|S|}\) associated with constraints (17) that were computed when solving the last relaxed CP (16)–(19) (i.e., without constraints (20)). Recall that, for each column index \(j\in S\), there are \(|T_j| -1\) constraints (17) in the CP, where \(T_j\) is the subset of tasks covered by \({\mathbf {A}}_j\). These constraints are associated with the first \(|T_j| -1\) tasks of \(T_j\). We denote by \({\hat{\pi ^q_j}}\), for all \(j\in S\) and \(q \in \{ 1, \ldots , |T_j| - 1\}\), the component of \({\hat{{\varvec{\pi }}}}\) associated with the task q covered by \({\mathbf {A}}_j\). Similarly, we denote by \(\alpha ^q_j\), for all \(j\in S\) and \(q \in \{ 1, \ldots , |T_j|\}\), the component of \({\varvec{\alpha }}\) associated with the task q covered by \({\mathbf {A}}_j\). To determine values for the \(\alpha ^q_j\) variables, we solve the following linear system of equations:

which can be decomposed by column index \(j\in S\). Conditions (21) ensure that the columns in the current solution of the SPP have a zero reduced cost and all compatible columns in the RP have, thus, a nonnegative reduced cost. Furthermore, conditions (22) ensure that the least value among all weighted reduced costs \({\bar{c}}_j / w_j = (c_j - {\varvec{\alpha }}^{\top } {\mathbf {A}}_j) / w_j\), \(j\in I_S\), of the incompatible columns considered in the CP is maximized. A dual solution \({\varvec{\alpha }}\in {\mathbb {R}}^m\) that satisfies conditions (21) is said to correspond to the current solution S.

Finally, observe that, if \(({\hat{{\varvec{\pi }}}}, {\hat{\sigma }}) \in {\mathbb {R}}^{m-|S|} \times {\mathbb {R}}\) is an optimal dual solution to the relaxed CP (16)–(19), then \(({\varvec{\alpha }}, {\hat{\sigma }}) \in {\mathbb {R}}^m \times {\mathbb {R}}\) forms an optimal solution to the relaxed CP (7)–(10) when \({\varvec{\alpha }}\) is computed by solving the equation system (21)–(22). In these solutions, \({\hat{\sigma }}\) is the dual value associated with constraints (9) and (18), which is equal to the optimal value of the relaxed CP.

Reprints and permissions