Daniel Cosmin Porumbel

François Clautiaux

Convergent Dual Bounds Using an Aggregation of Set-Covering Constraints for Capacitated Problems

Extended formulations are now widely used to solve hard combinatorial optimization problems. Such formulations have prohibitively-many variables and are generally solved via Column Generation (CG). CG algorithms are known to have frequent convergence issues, and, up to a sometimes large number of iterations, classical Lagrangian dual bounds may be weak. This paper is devoted to set-covering problems in which all elements to cover require a given resource consumption and all feasible configurations have to verify a resource constraint.

We propose an iterative aggregation method for determining convergent dual bounds using the extended formulation of such problems. The set of dual variables is partitioned into k groups and all variables in each group are artificially linked using the following groupwise restriction: the dual values in a group have to follow a linear function of their corresponding resource consumptions. This leads to a restricted model of smaller dimension, with only 2k dual variables. The method starts with one group (k = 1) and iteratively splits the groups. Our algorithm has three advantages: (i) it produces good dual bounds even for low k values, (ii) it reduces the number of dual variables, and (iii) it may reduce the time needed to solve sub-problems, in particular when dynamic programming is used.

We experimentally tested our approach on two variants of the cutting-stock problem: in many cases, the method produces near optimal dual bounds after a small number of iterations. Moreover the average computational effort to reach the optimum is reduced compared to a classical column generation algorithm.

Introduction

Column Generation (CG) is a well known technique for solving Linear Programs (LPs) with prohibitively-many decision variables. At each iteration, a CG method considers a Reduced Master Problem (RMP) with a limited number of variables. This RMP is solved to optimality so as to generate: (i) a primal bound for the CG optimum, and (ii) dual values that are provided as an input of a pricing sub-problem that computes a new variable of attractive reduced cost. By iteratively solving the sub-problem, columns are added to the RMP and the primal bounds (RMP optima) converge toward the original CG optimum. A dual bound can be determined at each iteration by computing the Lagrangian bound associated to the current dual multipliers.

A frequent issue of column generation methods concern their convergence and (lack of) stabilization. In numerous cases, a large number of iterations might be needed to obtain a useful dual bound. The last decades have seen a surge of interest in stabilization methods that aim at decreasing the number of iterations in CG [START_REF] Marsten | The BOXSTEP method for large-scale optimization[END_REF][START_REF] Du Merle | Stabilized column generation[END_REF][START_REF] Luebbecke | Selected topics in column generation[END_REF][START_REF] Carvalho | Using extra dual cuts to accelerate column generation[END_REF][START_REF] Oukil | Stabilized column generation for highly degenerate multiple-depot vehicle scheduling problems[END_REF][START_REF] Briant | Comparison of bundle and classical column generation[END_REF][START_REF] Clautiaux | New stabilization procedures for the cutting stock problem[END_REF].

Aggregation and disaggregation techniques have a long history in optimization [START_REF] Shetty | Solving large-scale linear programs by aggregation[END_REF]. Generally speaking, the goal is to transform programs with high degree of detail into more coarser programs of smaller size. For instance, one can aggregate close time instants [START_REF] Macedo | Solving the vehicle routing problem with time windows and multiple routes exactly using a pseudo-polynomial model[END_REF], nearby locations, related scheduling tasks, etc.; numerous examples and references are available in [START_REF] Litvinchev | Aggregation in large-scale optimization[END_REF]. In the remaining, we focus however on the specific context of column generation. Constraint aggregation is used in Figure 1: Our convergent method uses an "inner approximated dual polytope" that is iteratively refined (P 1 , P 2 , P 3 , . . . ) so as to reach the optimum of an initial polytope P with prohibitively many constraints. The optimum of any intermediate P k is a lower (dual) bound for the P optimum. column generation for stabilization reasons, to reduce the number of dual variables, or to produce smaller RMPs with less degeneracy [START_REF] Elhallaoui | Multi-phase dynamic constraint aggregation for set partitioning type problems[END_REF][START_REF] Elhallaoui | Dynamic aggregation of setpartitioning constraints in column generation[END_REF][START_REF] Benchimol | Stabilized dynamic constraint aggregation for solving set partitioning problems[END_REF].

We propose a new aggregation technique sharing certain similarities with the Dynamic Constraint Aggregation (DCS) method [START_REF] Elhallaoui | Multi-phase dynamic constraint aggregation for set partitioning type problems[END_REF][START_REF] Elhallaoui | Dynamic aggregation of setpartitioning constraints in column generation[END_REF], and, implicitly, with the Stabilized DCS (SDCS) method [START_REF] Benchimol | Stabilized dynamic constraint aggregation for solving set partitioning problems[END_REF] which combines DCS with a penalty-based stabilization technique [START_REF] Oukil | Stabilized column generation for highly degenerate multiple-depot vehicle scheduling problems[END_REF]. To construct new columns, DCS first disaggregates the dual values to obtain a dual solution expressed in the original dual space (see also Section 2.2). Then, the initial subproblem is solved and the resulting column is either added to the current aggregated RMP, or put aside if it is not compatible with an associated partition. In our method, we decide a priori how the dual variables will be disaggregated by adding cuts that link the dual variables in each group. We add to the RMP all columns that are computed by the aggregated subproblem and some of the resulting columns are only feasible in the aggregated model. This may lead to non-valid dual cuts with respect to the original dual polytope, i.e., excessively strong constraints. In other words, while DCS aims at reducing the number of iterations needed to converge while staying primal feasible, our method aims at producing useful iterative dual bounds by computing dual feasible solutions.

Such non-valid dual cuts generate, for each given number of groups k, an inner approximation P k of the dual polytope P. This approximated polytope P k has only two variables per group, and leads to an easier optimization problem when k is small. Once an optimal solution for P k is found, a valid dual bound is obtained for the original problem. The method converges toward the optimum by iteratively breaking groups until the optimum of the aggregated model corresponds to the optimum in the original dimension of P.

Fig. 1 presents an intuitive illustration of this convergent method. It constructs fast effective dual bounds in the first steps (coarse aggregation for small k) and performs incremental iterative calculations to progressively increase k. Compared to other aggregation methods, our approach has the advantage of providing valid dual bounds before fully converging, i.e., each P k optimum represents a valid dual bound for the original problem.

In this paper, we focus on set-covering master problems and pricing sub-problems modeled as Integer Linear Programs (ILPs) with a resource constraint. In these models, each element to cover requires a resource consumption. Our new dual groupwise-linearity restrictions work as follows: in each of the k groups, the dual values have to follow a linear function of corresponding resource consumptions. This allows us to reformulate the simplified problem with a new model whose variables are the slope and y-intercept of the linear function over each group. The motivation for using such aggregation comes from the fact that dual values are often correlated with the resource consumption at optimality.

The remainder is organized as follows. Section 2 presents the classical column generation method and a review of the literature. Section 3 is devoted to the convergent method: the first part ( §3.1- §3.3) describes the aggregation modelling for a fixed k and the second part presents the incremental calculations of the P 1 , P 2 , . . . sequence ( §3.4- §3.6). In Section 4, we discuss the theoretical benefits of using groupwise-linearity constraints instead of simpler equality constraints. Section 5 specifies the details of an application of the method to the cutting-stock problem. Section 6 presents our computational experiments, followed by conclusions in the last section.

2 Context, Motivation and Literature Review

Column Generation for Set-Covering Formulations with Resource Constraints in the Subproblem

We first introduce the set-covering model that is considered in this paper. Such models are often used to solve, among many others, cutting and packing [START_REF] Gilmore | A linear programming approach to the cutting stock problem[END_REF][START_REF] Gilmore | A linear programming approach to the cutting stock problempart II[END_REF], vehicle-routing [START_REF] Azi | An exact algorithm for a single-vehicle routing problem with time windows and multiple routes[END_REF], and employee scheduling [START_REF] Elhallaoui | Dynamic aggregation of setpartitioning constraints in column generation[END_REF][START_REF] Elhallaoui | Multi-phase dynamic constraint aggregation for set partitioning type problems[END_REF] problems.

The primal ILP is defined on a ground set of elements I = {1, . . . , n}. A column, also called configuration hereafter, of the master problem is an integer column vector a = [a 1 . . . a n ] ⊤ ∈ Z n + . The set A of all configurations is defined by all the solutions of a given sub-problem. The setcovering problem requires finding the minimum number of configurations that need to be selected so as to cover each element i ∈ [1..n] at least b i times (b i ∈ Z + ). To each configuration a ∈ A, we associate a variable λ a that indicates the number of times configuration a is selected. This leads to the following ILP: We are interested in the linear relaxation of this model, where the last constraint is replaced by λ a ≥ 0, ∀a ∈ A. The dual of this linear relaxation can be written using a vector y = [y 1 y 2 . . . y n ] ⊤ ∈ R n + of dual variables and a possibly exponential number of dual constraints. max b ⊤ y a ⊤ y ≤ 1, ∀a ∈ A y i ≥ 0, i ∈ 1, . . . , n P (2.2)

In the remainder, we will refer to this problem as the Dual Set-Covering Problem (DSCP) over polytope P: DSCP(P) = max{b ⊤ y : y ∈ P}. The optimum of this program, hereafter noted OPT(DSCP(P)), represents the main focus of all methods from this paper. It is also termed the CG optimum OPT CG , or the optimum over P on objective function b (we sometimes omit mentioning b, as we never use other functions).

Given the prohibitive large size of A, this constraint set is generated iteratively. From a dual point of view, the classical column generation can be seen as a cutting-plane algorithm: (i) reach the best dual solution in the current dual polytope defined by the limited set of constraints that have been generated so far, (ii) find a valid dual constraint that is violated by the current solution (report optimum if there is no such violated constraint), and (iii) refine the current formulation by adding a violated dual constraint and repeat from (i). Therefore, the column generation method constructs at each step an outer approximation of the dual polytope that is optimized.

The pricing sub-problem in CG requires finding a configuration of most negative reduced cost:

min a∈A 1 -y ⊤ a = 1 -max a∈A y ⊤ a (2.3)
In the following, we will forget the constant and seek a configuration a * that maximizes y ⊤ a * . Note that in this sub-problem, y represents input data; at each iteration, y is the dual optimum solution of the RMP.

A key point in column generation is the speed of the sub-problem solution method. Dynamic Programming (DP) algorithms provide an efficient tool for generating columns for many different families of problems (see among many others [START_REF] Azi | An exact algorithm for a single-vehicle routing problem with time windows and multiple routes[END_REF][START_REF] Gilmore | A linear programming approach to the cutting stock problem[END_REF]). We now formalize a general DP framework for problems with a single resource constraint. We consider a resource vector w = [w 1 w 2 . . . w n ] ⊤ ∈ R n + and a lower and upper limit of total resource consumption (C -, C + ∈ N). For a given dual solution y ∈ R n , the sub-problem requires selecting a i times each element i ∈ [1..n] (i.e., finding a configuration) so as to maximize the total profit y ⊤ a under the resource constraint:

C -≤ w ⊤ a ≤ C + .
The dynamic programming recursion computes a profit function P max that maps any state (c, i)

∈ [0, C + ] × [1.
.n] to the maximum profit P max (c, i) that can be obtained with a resource amount of c, only selecting elements of [1..i]. Generally speaking, this function can be determined using a recursion of the form below.

P max (c, i) = max P max (c -r • w i , i -1) + r • y i : r ∈ [0..b i ], r • w i ≤ c (2.4)
The initial state is P max (0, 0) = 0 and the final minimum reduced cost is attained in a state

(c * , n) such that P max (c * , n) = max c∈[C -,C+] P max (c, n).
If the given capacitated problem imposes additional conditions on feasible patterns, this formula is generalized by imposing additional constraints on the feasible subsets of elements. For example, in Bin Packing with conflicts [START_REF] Gendreau | Heuristics and lower bounds for the bin packing problem with conflicts[END_REF], some pairs of elements cannot be selected together.

The above recursion can be computed in O(n b C + ) time, where n b = n i=1 b i is the number of individualized elements. We use n b instead of n: the elements with demand multiplicities b i > 1 are not considered only once in the DP scheme, but b i times-observe r ∈ [0..b i ] in (2.4). For the Bin-Packing problem, n b = n and the complexity becomes O(nC + ). In other resource-constrained problems, such as the vector-packing problem [START_REF] Caprara | Lower bounds and algorithms for the 2-dimensional vector packing problem[END_REF], additional constraints can restrict the set of feasible patterns and/or increase the dimension of the state space.

Although we describe here the case with one resource, it can also be applied when several resources are considered. One only needs multiple dimensions to index the states.

Aggregation in Column Generation and Extended Formulations

A way to cope with prohibitively large Mixed-Integer Programs (MIPs) arising in extended formulations is to approximate them. This can be done by restricting the MIP to a subset of variable and/or constraints. This leads to a primal inner approximation that can be refined iteratively (typically by column-and-row generation, see [START_REF] Sadykov | Column generation for extended formulations[END_REF] for a generic view of those methods).

Another way of obtaining a tractable model is to apply aggregation to the constraints or to the variables. This is done statically by [START_REF] Vyve | Approximate extended formulations[END_REF], who define a smaller model whose size depends on a given parameter. An interesting conclusion of [START_REF] Vyve | Approximate extended formulations[END_REF] is that small values of this parameter are often sufficient to obtain excellent bounds. However this method is static and does not converge toward the optimum of the initial model.

In the context of column generation methods, similar aggregations can be obtained by adding equality constraints on dual variables of a same cluster. This has the advantage to reduce the size of the pricing sub-problem (since all dual variables of the cluster can be exchanged), and also to stabilize the column generation process. This strategy is useful if many dual variables are "equivalent" and thus many of them can have similar values at optimality. Dual optimal and deep dual optimal inequalities [START_REF] Amor | Dual-optimal inequalities for stabilized column generation[END_REF] are examples of dual constraints that do not separate all optimal dual solutions.

Dynamic constraint aggregation has been introduced in the context of column generation in [START_REF] Elhallaoui | Dynamic aggregation of setpartitioning constraints in column generation[END_REF][START_REF] Elhallaoui | Multi-phase dynamic constraint aggregation for set partitioning type problems[END_REF][START_REF] Amor | Dual-optimal inequalities for stabilized column generation[END_REF]. In the aggregated model, each constraint represents a cluster of original partitionning constraints. From a dual space perspective, a dual variable of the aggregated model represents a cluster of original dual variables, and therefore its value is equal to the sum of the actual dual values in the cluster. When the pricing subproblem is called, the dual variables are disaggregated and the initial subproblem is solved. The column produced by the subproblem is added to the RMP if it is compatible with the current aggregation, or put aside otherwise. These methods ensure the exactness of the solution by iteratively breaking up the clusters until an optimal feasible solution is obtained.

In the above approach, a compatible column corresponds to a configuration that either covers all elements of a given group or none of them. The elements of a given group are not distinguished in the aggregated model, i.e., they are in the same equivalence class and they are all associated to a unique representative partitioning constraint. For some specific problems, such as the cuttingstock problem, we are given items with different weights (resource consumptions); if the weights of two items are different, they might have significantly different dual values at optimality, and so, they can hardly be considered equivalent. In such situations, equivalence-based approaches might require too many (small) clusters to obtain good approximations. An extreme case is obtained when each cluster must contain exactly one element, even if the structure of the solution is simple (for example the dual values can be obtained by applying a linear function to the item sizes).

In certain cases, the dual values may often follow a piecewise-linear function of the associated resource consumptions. This type of solution can be obtained by a new type of aggregation, where dual variables in the same cluster have values that follow a linear function of this resource. This is what we intend to do in the remainder of this paper. The basic underlying idea is that elements requiring higher resource consumptions arise in less valid configurations, and so, they can be associated to larger dual values; for cutting-stock, the dual solution vector is always nondecreasing (see for example [START_REF] Carvalho | Using extra dual cuts to accelerate column generation[END_REF]).

Iterative Inner Dual Approximation (2IDA)

This section describes the iterative aggregation method, starting from a general overview and gradually addressing all technical details. From now on, we will focus on the the dual polytope P, i.e we consider CG as a cutting-plane algorithm for optimizing (2.2).

Full Description of the Aggregated Dual Polytope

Let G k = {I 1 , I 2 , . . . I k } be a partition of I = {1, 2 . . . n} into k groups, each group having n j = |I j | elements. Without loss of generality, we consider that the elements in each group are sorted (indexed) in increasing resource consumption. Given elements I j of a group j ∈ [1..k], let y j , w j , b j and a j denote the n j -dimensional column vectors related to dual variables, resource consumptions, demands and, respectively, coefficients of some configuration a ∈ A. Variables in the dual vector are ordered in such a way that y is the concatenation of k group components:

y ⊤ = [y 1 y 2 . . . y n ] = [(y 1 ) ⊤ . . . (y k ) ⊤ ].
We artificially impose the following groupwise-linearity restrictions: given any group j ∈ [1.

.k], the values of the dual variables y j 1 , y j 2 , . . . y j nj have to follow a linear function of the resource consumption of the corresponding element. This function depends on a slope α j and a y-intercept β j . More precisely, the i th component of y j can be written as y j i = w j i α j + β j . By simply adding these constraints to (2.2), we obtain the inner dual polytope:

max b ⊤ y a ⊤ y ≤ 1, ∀a ∈ A y j i = w i α j + β j , ∀j ∈ [1..k], i ∈ I j y j i ≥ 0, ∀j ∈ [1..k], i ∈ I j α j , β j ∈ R, ∀j ∈ [1..k] (3.1)
Remark 1. By extracting the vector y from any valid solution of above (3.1), we obtain a valid solution of P in (2.2). The opposite is not true in general: only vectors y of P with the suitable groupwise linear structure can be lifted to valid solutions of (3.1).

Let us re-write (3.1) using only the new decision variables α j and β j . We first re-write the objective function. For each group j, y j can be written as a linear combination of w j and 1 j (vector [1 1 . . . 1] ⊤ with n j elements): y j = w j α j + 1 j β j . Observe that

b ⊤ y = k j=1 b j ⊤ y j = k j=1 b j ⊤ w j α j + 1 j β j = k j=1 b j ⊤ w j α j + b j ⊤ 1 j β j (3.2)
The first set of constraints from (3.1) can be written k j=1 a j ⊤ y j ≤ 1, a ∈ A. For each j, we have:

a j ⊤ y j = a j ⊤ w j α j + 1 j β j = a j ⊤ w j α j + a j ⊤ 1 j β j (3.3)
We are now ready to express model (3.1) with variables α j and β j only. To simplify the notation, we introduce the following definitions.

Definition 1. Given a configuration a ∈ A and a group j, we define:

-c j a = a j ⊤ w j : total resource consumption of the elements of I j selected in a. Remark this is the coefficient of variable α j in (3.3); -N j a = a j ⊤ 1 j : total number of elements of I j selected in a. This is the coefficient of variable β j in (3.3); -w j min and w j max : the lowest and respectively highest resource consumption of an element of group j.

By simply substituting (3.2)-(3.3) in model (3.1), after handling the non-negativity constraints, we obtain an equivalent model in the space R 2k :

max k j=1 b j ⊤ w j α j + b j ⊤ 1 j β j k j=1 c j a α j + N j a β j ≤ 1, ∀a ∈ A w j min α j + β j ≥ 0 ∀j ∈ [1..k] w j max α j + β j ≥ 0 ∀j ∈ [1..k] α j , β j ∈ R, ∀j ∈ [1..k]        P k (3.

4)

Proposition 1. There is a bijection between the set of feasible values of [y 1 y 2 . . . y n ] in the model

(3.1) in R n + × R 2k and the feasible solutions [α 1 . . . α k β 1 . . . β k ] ⊤ of model (3.4) in R 2k .
Proof. Recall that (3.1) has by construction the constraint y j i = w i α j + β j . The bijection simply maps any solution [y 1 . . . (3.4). Recall (3.4) is simply obtained by replacing all y j i terms of (3.1) with y j i = w i α j + β j , using (3.3) and (3.2). The non-negativity constraints of (3.1) are simplified in (3.4), i.e., y j i = w j i α j + β j ≥ 0, ∀j ∈ [1.

y n , α 1 . . . α k , β 1 . . . β k ] ⊤ ∈ R n + × R 2k that is valid in (3.1) to a solution [α 1 . . . α k , β 1 . . . β k ] ⊤ ∈ R 2k that is valid in
.k], i ∈ I j is satisfied if y j min = w j min α j + β j ≥ 0 and y j max = w j max α j + β j ≥ 0: the value w j i α j is always between w j min α j and w j max α j . The new formulation (3.4) has 2k real variables that can be either positive or negative. In this basic version, only the number of variables is decreased. This reduction can be useful for stabilization reasons in CG: using less dual variables, one reduces the possibility of oscillation. Also, less dual variables corresponds to less constraints in the corresponding Reduced Master Problems (RMPs), which can accelerate the LP solver for RMPs.

However, a drawback of (3.4) consists of the fact that it has the same (prohibitively large) number of constraints as (2.2). We address this point in the next section, introducing a new model (3.5) that removes many constraints that become redundant in (3.4).

Column Generation for the Agregated Set-Covering Problem

When one imposes linearity restrictions on the dual variables, many configurations of A are not useful in (3.4) anymore. This will be shown in Section 3.2.1 below, leading to a new model (3.5) that reduces the size of the initial set-covering LP in terms of both columns and rows. This latter model is the one which is actually optimized by our method; more computational and optimization discussion follow in Section 3.2.2.

Removing Redundant Constraints in Model (3.4)

Model (3.4) does not explicitly indicate the non-zero elements in a configuration a ∈ A. The main constraints of (3.4) are built only by considering the number of elements of each group N j a and the total resource consumption c j a . Therefore, to avoid redundancies, we no longer express configurations as vectors in R n + , but as aggregated R 2k vectors of the form a

= [c 1 a N 1 a , c 2 a N 2 a , . . . , c k a N k a ]
. We say that the aggregated configuration a can be lifted to a configuration a ∈ A if:

1. C -≤ k j=1 c j a , ≤ C + , i.e.
, the total consumption is feasible;

2. there exist a ∈ A such that c j a = a j ⊤ w j and N j a = a j ⊤ 1 j , for all j ∈ [1..k].

There might exist more than one disaggregated configuration a ∈ A that corresponds to an aggregated configuration a, i.e., more configurations a might have the same values c j a , N j a , ∀j ∈ [1..k]. Therefore, many constraints of P k in (3.4) are redundant. Definition 2. Given group j ∈ [1..j], the set R j of feasible resource consumptions is: R j = c j ∈ [0..C + ] : ∃a ∈ A such that c j = c j a Definition 3. Given a feasible consumption c j ∈ R j , we define the minimum and maximum cardinality coefficients of c j for group j:

-N + (j, c j ): the maximum number of I j elements in a valid configuration a ∈ A in which the total resource consumption of the I j elements is c j ; -N -(j, c j ): the minimum number of I j elements in a valid configuration in which the total resource consumption of the I j elements is c j . Definition 4. Let A k be the set of dominant (non-redundant) configurations, defined by:

A k = a ∈ A, such that ∀j ∈ [1..j], N j a ∈ N + (j, c j a ), N -(j, c j a )
By replacing A with A k in model (3.4), we obtain a new model:

max k j=1 b j ⊤ w j α j + b j ⊤ 1 j β j k j=1 c j a α j + N j a β j ≤ 1, ∀a ∈ A k w j min α j + β j ≥ 0 ∀j ∈ [1..k] w j max α j + β j ≥ 0 ∀j ∈ [1..k] α j , β j ∈ R ∀j ∈ [1..k]        P k (3.5)
Obviously, if two configurations a, a ′ ∈ A lead to the same values of N j a and c j a , only one constraint has to be explicitly considered. The new optimization problem, will be hereafter referred to as the Dual Set-Covering Problem (DSCP) over polytope

P k : DSCP(P k ) = max{ k j=1 ((b j ) ⊤ w j )α j + ((b j ) ⊤ 1 j )β j : [α 1 . . . α k , β 1 . . . β k ] ⊤ ∈ P k }.
We will show below that (3.5) determines the same polytope P k as (3.4), and so, any solution of (3.5) can be lifted to a solution of the initial LP (2.2).

Proposition 2. Any solution of (3.5) can be lifted to a solution that is feasible for the initial set-covering model (2.2).

Proof. Remark 1 (p. 5) states that any solution of P k in (3.1) can be projected into a feasible solution of the initial dual polytope P in (2.2). Proposition 1 shows that P k is modelled equivalently by (3.1) and (3.4). It is enough to show that (3.4) is equivalent to (3.5).

Using the notation from Definitions 1 and 3, we observe that, given any valid a ∈ A, for each j ∈ [1..k], the values of variables α j and β j have to respect one of the two inequalities below:

1. c j a α j + N j a β j ≤ c j a α j + N + (j, c j a )β j , if β j ≥ 0; 2. c j a α j + N j a β j ≤ c j a α j + N -(j, c j a )β j , if β j < 0.
Therefore, the constraints obtained from configurations outside A k are dominated by constraints associated to configurations a ∈ A k , i.e., with N j a = N + (j, c j a ) or N j a = N -(j, c j a ). All configurations outside A k are dominated, and so, (3.4) is equivalent to (3.5).

To solve DSCP(P k ), we need to determine all coefficients N + (j, c j ) or N -(j, c j ), ∀c j ∈ R j , j ∈ [1..k]. This is performed in a coefficient calculation preprocessing stage which is executed only once before launching column generation. The goal is to find the maximum and minimum number of elements that can be selected from I j so as to consume a total resource amount of c j , ∀c j ∈ R j , j ∈ [1..k]. This task is similar to the pricing sub-problem (2.3) for the non-aggregated problem (Section 2.1). If the dynamic programming scheme for (2.4) is fast, so is the dynamic programming for this coefficient calculation. Indeed, it is enough to replace y i with 1 in (2.4) and to apply an analogous dynamic programming twice, i.e., once with a minimization and once with a maximization objective. Such a dynamic programming scheme generates by default a state for each feasible c j value, i.e., the set R j . Considering all groups j ∈ [1..k] together, this coefficient calculation stage requires the same complexity as one pricing sub-problem for the initial model (2.2).

Remark 2. The number of constraints in (3.5) is exponentially large in k, but not in n. Computing all constraints of P k would thus require an asymptotic running time that is exponentially large in k and polynomial in C + (to determine all N + (j, c j ) and N -(j, c j )). If we consider that k and C are bounded constants, it is possible to achieve the full construction of P k in polynomial time. This would lead to a pseudo-polynomial dual bound for the initial combinatorial optimization problem. Although such direct P k constructions are used in this paper only for k = 1, the main idea of a direct construction represents a promising subject for further work.

An aggregated Pricing Subproblem

The model (3.5) is optimized by column generation for all values of k > 1. A computational difficulty is related to the pricing algorithm. In a straightforward approach, one could disaggregate the variables to obtain the dual solution y in the original space in R n + ; then, the pricing algorithm for the non-aggregated model (Section 2.1) can be used. This has the disadvantage of considering a solution space with many symmetries and redundancies, since all y j i with i ∈ [1..n j ] in a group j ∈ [1..k] correspond to a single pair of dual variables α j and β j .

We now show how the pricing can be solved without disaggregation. Iteratively, for each current solution [α 1 . . . α k , β 1 . . . β k ] ⊤ ∈ P k , the column generation algorithm needs to solve an aggregated version of the sub-problem (2.3): find the aggregated configuration a

= [c 1 a N 1 a , c 2 a N 2 a , . . . , c k a N k a ] that maximizes the profit value k j=1 α j c j a + β j N (j, c j a )
. As for the non-aggregated pricing (2.3), a column of negative reduced cost is identified whenever this profit is larger than 1. Formally, the pricing problem is the following:

max k j=1 α j c j + β j N (j, c j ) s.t. C -≤ j j=1 c j ≤ C + N (j, c j ) = N + (j, c j ) if β j ≥ 0 N -(j, c j ) if β j < 0 ∀j ∈ [1..k] c j ∈ R j , ∀j ∈ [1..k] , (3.6)
where R j is the set of feasible resource consumptions (Definition 2). This set is determined during the preprocessing stage for coefficient calculation from Section 3.2.1, i.e., R j is implicitly determined at the same time with the values N + (j, c j ) and N -(j, c j ). Recall that the dynamic programming approach proposed for this coefficient calculation requires the same asymptotic running time as one pricing algorithm for the original problem (2.2), for all j ∈ [1.

.k].

The variables c j ∈ R j , ∀j ∈ [1.

.j] are actually the only decision variables in (3.6). They are sufficient to represent a solution for this pricing sub-problem, since any N (j, c j ) term can be deduced from c j and β j . All β j represent input data for the sub-problem, and so, we choose from the beginning to use either N (j, c j ) = N + (j, c j ) or N (j, c j ) = N -(j, c j ), depending on the sign of β j . Let us denote p(j, c) = α j c + β j N (j, c) the potential profit that can be obtained with a resource amount of c for group j. The pricing problem (3.6) is reformulated:

max k j=1 c∈R j p(j, c)x jc s.t. C -≤ j j=1 c∈R j c j x jc ≤ C + , c∈R j x jc = 1, ∀j ∈ [1..k] x jc ∈ {0, 1}, ∀j ∈ [1..k], c ∈ R j (3.7)
where x jc is a new decision variable for group j, i.e., it is equal to 1 if group j has a total resource consumption of c, and 0 otherwise.

The resulting aggregated pricing is a multiple-choice variant of the non-aggregated pricing from Section 2.1. For instance, if the non-aggregated problem is the knapsack problem, the aggregated one is the multiple-choice knapsack problem. Generally speaking, the non-aggregated DP approach from Section 2.1 can be extended to an aggregated Dynamic Program (DP) as follows. We associated a DP state to each total consumption value in [1..C + ], and, for each level j ∈ [1.

.k], one can scan all feasible choices (total consumptions values c j ) to generate (or update) other states. Furthermore, in certain cases, the aggregated DP only needs to perform these computations on two levels associated to the new sub-groups j 1 , j 2 ∈ [1..k] (see Section 3.6), which is fast in practice. However, in this aggregated pricing, the number n b = n i=1 b i of individualized elements is no longer a factor in the asymptotic or practical running time, once the initial preprocessing (Section 3.2.1) is performed.

3.3

The Iterative Algorithm and the Incremental P k Optimization Section 3.2 described an aggregation-based column generation method that optimizes (3.5) for a given fixed partition G k = {I 1 , I 2 , . . . I k }. The resulting value, hereafter noted lb G k (or simply lb k when the exact structure of G k is not essential) represents a lower bound for the sought OPT CG . Algorithm 1 below presents the general steps of our Iterative Inner Dual Approximation (2IDA) method. The main idea is to iteratively break the groups into smaller subgroups and incrementally refine the description of P k until a stopping condition is met.

Algorithm 1: Iterative Inner Dual Approximation k ← 1, G 1 ← {I} 1 P k ← detP k Coefs(G k ) // preprocessing coefficient calculation (Section 3.2.1) 2 repeat 3 lb G k ←aggregCGoptim(P k ) // Section 3.2 (or full P k construction for k = 1) 4 ub G k ← upBound(P k , lb G k ) // optional upper bound (Section 3.5) of OPT CG 5 if ⌈lb G k ⌉ < ⌈ub G k ⌉ then 6 G k+1 ← grpSplit(G k )
// split a group and refine partition (Section 3.4)

7 P k+1 ← detP k Coefs(G k+1 -G k ) //determine
9 k ← k + 1 10 until (⌈lb G k ⌉ = ⌈ub G k ⌉) or a stopping condition is reached ; 11
The partition G k of the dual variables is initialized with one group equal to I. Right after that (in Line 2), 2IDA applies the preprocessing stage from Section 3.2.1 to determine the coefficients (e.g., N -(j 1 , c j1 ), N + (j 1 , c j1 ), etc.) required to describe all aggregated polytopes. Routine aggregCGoptim(P k ) uses partition G k to restrict and project the dual polytope P to dimension R 2k and describe polytope P k (3.5). This polytope has fewer constraints and fewer variables than P. For k > 1, aggregCGoptim(P k ) relies on a column generation scheme that iteratively solves the aggregated pricing sub-problem from Section 3.2.2.

After determining lb G k , Algorithm 1 can continue by calculating an optional upper bound ub G k (routine upBound(P k , lb G k ), Line 5). For this, we propose (see Section 3.5) a method that performs disaggregated CG steps, using a polytope P u that contains much fewer constraints than P. While upper bounds are not mandatory in 2IDA, they can be useful for "tail cutting" (i.e., to stop 2IDA if ⌈lb G k ⌉ = ⌈ub G k ⌉). Furthermore, the specific upper bound described in Section 3.5 can also be more effective in guiding group split decisions (see below).

Section 3.4 presents several methods for taking group split decisions (call grpSplit(G k ) in Line 7). The goal is to determine a group j * ∈ [1..k] to be split into sub-groups j 1 and j 2 , so as to progress from k to k + 1. This decision can be taken using information related to the upper bound (if available), problem-specific data, or general ad-hoc indicators (group sizes, resource consumption spreads, etc.). After calling grpSplit(G k ), 2IDA applies (in Line 8) the coefficient calculation preprocessing stage (Section 3.2.1), but only restricted to the new sub-groups j 1 and j 2 .

The configurations (constraints) generated for a given k are not completely removed when 2IDA passes to k + 1. The general goal of 2IDA is to perform an incremental construction of P 1 , P 2 , P 3 . . . . Instead of constructing each P k+1 from scratch, 2IDA picks up constraints generated when optimizing over P k and lifts them in the new dimension of P k+1 . In fact, routine inheritCnstr(P k ,G k+1 ) make use of P k constraints to warm-start the optimization over DSCP(P k+1 )-see Section 3.6 for a description in greater detail of this "inheritance" process.

Split Decisions

2IDA converges towards the optimum of the initial problem for any group split decision. However, from a practical perspective, the split decision is a crucial part for the efficiency of the method. In many optimization settings, a fast convergence is obtained by making at each step a decision that leads to the best local improvement. However, the split decision that leads to the best improvement for a given step k can only be known after optimizing over P k+1 . Trying all possible split decisions and optimizing all possible DSCP(P k+1 ) would definitely require prohibitive calculation time. We present below techniques for making good splitting decisions using a limited computing time.

Splitting Groups Based on a Binary Decomposition

The most general splitting approach only works with regular interval lengths. To determine the group to be split, a reasonable choice consists of selecting a group j * with maximum resource consumption spread, i.e., such that w j * max -w j * min = max w j max -w j min . This group j * is split into sub-groups j 1 and j 2 using a split point c * ∈ [w j min ..w j max ] that separates the elements of I j * as follows:

I j1 = {i ∈ I j * : w i ≤ c * } with n j1 = |I j1 | and I j2 = {i ∈ I j * : w i > c * } with n j2 = |I j2 |.
To maintain the regularity of the intervals sizes, this split point is simply determined via c * = t 2 z C + , where z is the minimum integer that can lead to n j1 , n j2 ≥ 1 for some t ∈ Z + ; the best value of t is chosen afterwords so as to minimize |n j1 -n j2 |. This basic method is general, but does not take into account any information from the current step of the optimization.

Group-Splitting Following an Outside Reference Solution

The split method from this section aims at removing some P k constraints that block the optimum of DSCP(P k ) in (3.5) from advance towards an outside reference solution of better objective value. Definition 5. A solution y u ∈ R n + is an outside reference solution for P if b ⊤ y u ≥ OPT CG , i.e., y u is associated to an upper bound of OPT CG . This solution is either the optimal solution of (2.2), or belongs to the exterior of P.

We focus on the differences between the solution that realizes OP T (DSCP(P k )) and a given outside reference solution y u . In the following, we denote by y * k the optimal solution of DSCP(P k ) lifted in the original space R n + . The goal of this split method is to remove the constraints of P k that separate y * k from y u , and, in the first place, those which are verified to equality by y * k . In other words, we seek the constraints that should be removed to allow the solution y * k+1 of the next iteration to be closer to y u than y * k . We first introduce the notion of open direction. Definition 6. Let y be a feasible solution of P k and y u an outside reference solution. The direction y → y u is called "open" if there exists a small enough ǫ > 0 such that y ǫ = y + ǫ(y uy) ∈ P.

In concrete terms, there is a feasible solution y ǫ different from y that can be obtained by convex combination of y and y u .

We always use outside reference solutions y u generated from the current y * k by the upper bounding method from Section 3.5. We will prove (Proposition 3) that this only leads to reference solutions y u such that the direction y * k → y u is open (unless if b ⊤ y * k = OPT CG , in which case b ⊤ y u = OPT CG and 2IDA terminates). Let us focus on the constraints of P k that separate y * k from y ǫ = y * k + ǫ(y uy) for a small enough ǫ. Since y ǫ belongs to P, y * k and y ǫ are separated by a set of P k inequalities that do not belong to P. Such inequalities can only arise by combining valid P constraints with group-wise linearity restrictions associated to P k . We will have more to say about this process in Section 4, but for the moment is is enough to insist on the following. Dropping these non-valid restrictions could drive 2IDA towards y u , improving the objective value by at least some ǫb ⊤ (y uy).

Therefore, a split decision should select a group j * that is associated to groupwise linear restrictions separating y * k from y u . Since non-valid P k constraints spring from enforcing the dual values in each group to follow a linear function, we consider groups whose dual values in y u do not follow a linear function. Based on such information, we determine j * based on evaluating how far y u is from a linear function over each j ∈ [1..k]. An example of precise evaluation formula is provided in Section 5.3, i.e., function h 3 in the cutting stock application.

Such scores can be eventually combined with other indicators based on the resource consumption spread of the groups (as in Section 3.4.1). This allows one to assign a global improvement probability score for each possible split decision. This approach does not always lead to the direction of highest improvement, but it does use information on the current primal and dual bounds to guide the search. Moreover it does not have a high computational cost when such the outside reference solution y u is provided by the upper bounding approach of Section 3.5.

Upper Bounding and Open Directions based on the P k Description

The method proposed in this section takes y * k (optimal solution of DSCP(P k ) expressed in dimension R n + ) as input and returns an upper bound b ⊤ y u associated to an outside reference solution y u such that:

1. if b ⊤ y * k < OPT CG , the direction y * k → y u is open 2. when b ⊤ y * k = OPT CG , then b ⊤ y u = OPT CG .
The first property is useful for guiding the group splitting heuristic from Section 3.4.2 above. It means that P contains a linear combination of y * k and y u that is better than y * k . The second property allows to stop the 2IDA iterative process as soon as b ⊤ y * k = OPT CG .

3.5.1 Theoretical Modelling: the Polytope P u of y * k -Tight Constraints Given the optimal solution y * k ∈ R n + of OPT(DSCP(P k )) at current step k, we determine y u by optimizing the dual set-covering problem associated to polytope P u ⊃ P defined by:

P u = y ∈ R n + : a ⊤ y ≤ 1, ∀a ∈ A such that a ⊤ y * k = 1 (3.8)
We say that P u is the polytope of the y * k -tight constraints of P. It is simply obtained from P by only keeping the constraints that are satisfied to equality by y * k . The following proposition shows that P u ⊃ P, and so, OPT(P u ) ≥ OPT CG , i.e., y u is an outside reference solution of P that yields the upper bound value b ⊤ y u . Proposition 3. A solution y u such that b ⊤ y u = max{b ⊤ y : y ∈ P u } = OPT(DSCP(P u )), where P u is defined by (3.8), satisfies the following properties:

(1) upper bounding: b ⊤ y u ≥ OPT CG ;

(2) open direction: if y * k is not an optimal solution of DSCP(P), then b ⊤ y * k < b ⊤ y u and y * k → y u is an open direction;

(3) optimality proving: if b ⊤ y * k = OPT CG , then y u also satisfies b ⊤ y u = OPT CG . Proof. The first property comes directly from the P u definition (3.8). Since P u is constructed from a subset of the constraints of P, we directly have P u ⊃ P. By optimizing in direction b over P u , we obtain a solution y u that dominates the optimum over P in direction b, i.e., b ⊤ y u ≥ OPT CG .

We now prove [START_REF] Ben Amor | Cutting stock problems[END_REF]. Since y * k is not an optimal solution, we obtain b ⊤ y u ≥ OPT CG > b ⊤ y * k . Suppose that there is no feasible convex combination of y * k and y u besides y * k . In this case, y * k belongs to a facet of P, i.e., it cannot lie in the interior of P. This face would yield a y * k -tight constraint that would cut off all convex combination of y * k and y u , except y * k , but including y u . This is impossible because, by construction, y u satisfies all y * k -tight constraints. Therefore, the segment joining y * k and y u does contain more solutions that belong to P, and so, y * k → y u is an open direction.

The last property (3) comes from the construction of P u . If y * k is optimal in P, its basic constraints need to be y * k -tight, and so, they belong to P u . These constraints are sufficient to ensure the optimality of y u in P.

Practical Aspects: Constructing P u

We solve DSCP(P u ) by column generation. DSCP(P u ) has exactly the same structure as DSCP(P) in (2.2), but the set of possible constraints is restricted to those that are y * k -tight. The pricing sub-problem requires finding a configuration a ∈ A with a ⊤ y * k = 1 (y * k -tightness) that maximizes y ⊤ a for the current dual solution y. This pricing sub-problem can be solved exactly as the initial pricing sub-problem (Section 2.1), by replacing max a∈A y ⊤ a in (2.3) with a hierarchical objective function max a∈A (M y * k + y) ⊤ a, where M is a sufficiently large value. In practice, this behavior is obtained using a lexicographic maximization objective, i.e., first maximize a ⊤ y * k , and, subject to this, maximize a ⊤ y. The original dynamic programming scheme (Section 2.1) is not modified and the number of states is the same. Instead of computing a unique profit function P max in the recursion (2.4), the lexicographic dynamic program computes an ordered pair of profit functions (P * max , P u max ). For each state (c, i)

∈ [0, C + ]×[1.
.n], P * max (c, i) and P u max (c, i) represent the maximum of (y * k ) ⊤ a and, respectively, of y ⊤ a over all (partial) configurations a in state (c, i), i.e., consuming a resource amount of c only using elements of [1..i]. The recursion (2.4) evolves to a lexicographic maximization formula.

Proposition 4. Given a dual solution y, the above lexicographic dynamic program (optimizing (y * k ⊤ a, y ⊤ a) lexicographically) determines a configuration a ∈ A that maximizes y ⊤ a subject to (y * k ) ⊤ a = 1. Proof. The evolution of the lexicographic dynamic program can be interpreted as follows: at each recursion step, the lexicographic maximization formula actually maximizes (y * k ) ⊤ a, breaking ties using the value of y ⊤ a. As such, the y objective is only used for tie breaking, and so, the lexicographic dynamic program does return a configuration a * that maximizes the y * k objective.

As such, (y * k ) ⊤ a * achieves the maximum value of 1, because at least some "≤ 1" constraints in (3.5) need to be saturated by y * k (otherwise, no other (3.5) constraint can render infeasible a higher value solution such as y * k + ǫ1 n , with a sufficiently small ǫ).

Consequently, the complexity of the DP pricing for DSCP(P u ) is the same as the complexity of the DP pricing for DSCP(P). However, since P u is defined by a far smaller number of constraints than P, the proposed upper bounding method is generally faster.

P k+1 Optimization: Using P k to Warm-Start the Column Generation

After solving DSCP(P k ) to optimality at iteration k, 2IDA splits a group to obtain a new partition G k+1 with k + 1 groups. A new program (3.5), associated with a new polytope P k+1 has to be optimized. One could directly apply from scratch the column generation method from Section 3.2.2 after each k incrementation, but this may globally require a too large computational cost. Therefore, we use the columns already generated at iteration k to warm-start the column genphase at the next iteration k + 1. The main goal is to perform an incremental iterative description of the polytopes P 1 , P 2 , . . . . Let j * be the index of the group that is split after iteration k, and j 1 , j 2 be the two new (sub-)group indices. Right after taking the split decision, 2IDA (see Algorithm 1) starts building P k+1 . It first determines the new P k+1 coefficients (N -(j 1 , c j1 ), N + (j 1 , c j1 ), etc.) by calling detP k Coefs (Line 8). This routine applies on groups j 1 and j 2 the coefficient calculation preprocessing stage from Section 3.2.1. In fact, the coefficients of j 1 can be directly recovered from the program that was run on j * , supposing that the dynamic program for I j * considered the I j1 elements before I j2 .

Right after splitting j * , the next step of 2IDA (Line 9) generates some P k+1 constraints by lifting existing P k constraints from dimension R 2k to dimension R 2k+2 . The optimal solution of DSCP(P k ) is also lifted to construct an initial feasible solution for P k+1 , i.e., 2IDA keeps this solution unchanged over groups of [1..k] \ {j * }, and applies α j1 = α j2 = α j * and β j1 = β j2 = β j * for j 1 and j 2 . Starting from this P k+1 solution, 2IDA first applies column generation to optimize over an intermediate polytope (noted P ′ k+1 , see below) containing inherited constraints only. If the optimal solution of DSCP(P k ) remains optimal for DSCP(P ′ k+1 ), there is no need to continue optimizing over P k+1 . Otherwise, when no more inherited constraints lead to negative reduced cost columns, 2IDA proceeds to the full CG from Section 3.2.2, i.e., when Algorithm 1 goes on to Line 4.

More formally, we note A base k the set of configurations associated to constraints of P k that are either related to basic primal variables, or non-basic primal variables of reduced cost 0 (with regards to y * k ). We consider only the constraints that have been explicitly generated when optimizing P k in (3.5).

Our approach is equivalent to first constructing a dual polytope P ′ k+1 focusing on a set of inherited constraints:

A ′ k+1 = a ∈ A k+1 : ∃a ′ ∈ A base k , such that c j a = c j a ′ and N j a = N j a ′ ∀j ∈ [1.
.k] \ {j * } , where A k+1 is the set of (non-dominated) P k+1 constraints in (3.5).

Since A ′ k+1 ⊂ A k+1 , we have OPT(DSCP(P ′ k+1 ))) ≥ OPT(DSCP(P k+1 ))). Except for the split group j * , all coefficients of existing A base k configurations can be transmitted to recover more rapidly some A ′ k+1 configurations. Each time the pricing sub-problem is called, we seek and add negative reduced cost columns in two phases: (i) try to find columns associated to configurations in A ′ k+1 and (ii) search column using the original method from Section 3.2.2. For each configuration in A base k , step (i) seeks a configuration of negative reduced cost that uses the same coefficients (c j , N j ) for any group j = j * and new coefficients for j 1 and j 2 . This is solved using the multi-choice dynamic program from Section 3.2.2, but with only two decision levels j 1 and j 2 . This latter dynamic programming is significantly faster than the original one from Section 3.2.2. Besides only needing two levels instead of k, the maximum resource amount to consider is also reduced from C + to C + -j∈[1..k]-{j * } c j (the space occupied by inherited values). In certain cases, OPT(DSCP(P ′ k+1 )) = OPT(DSCP(P k )) and this is enough to conclude OPT(DSCP(P k+1 )) = OPT(DSCP(P k )), because the 2IDA bounds can not decrease by incrementing k. In practical terms, this fully justifies using the aggregated pricing approach from Section 3.2.2: when using this inheritance with two levels, it can sometimes incrementally compute lb k+1 from lb k in virtually no-time. The full algorithmic template is presented in Algorithm 2: Phase 1 roughly corresponds to Line 9 (in Algorithm 1) and Phase 2 to Line 4 (with incremented k). -Keep the values α j and β j for all j = j * // no change outside split group j * -α j1 ← α j , α j2 ← α j -β j1 ← β j , β j2 ← β j repeat for a ∈ A base k do -solve the aggregated multiple-choice pricing variant (Section 3.2.2) with 2 elements (j 1 and j 2 ) and capacity limits

C --j∈[1..k]-{j * } c j a and C + -j∈[1..k]-{j * } c j

Comparing Linear and Equality Aggregations

We now investigate how group-wise linear restrictions and valid P constraints produce together the "inner polytope" P k ⊂ P. We motivate on the use of linear aggregation with regard to simpler equality aggregations (e.g., y i = y j , for aggregated i, j ∈ [1..n]). Equality aggregations are mostly useful when one can assume that elements of similar weights will be associated to the same dual value in an optimal solution. This is sometimes true, for instance, in certain instances of the cutting-stock problem, if the weights of the aggregated items are similar. More generally, it is reasonable to perform equality aggregations on rather similar entities in the model, such as, related scheduling tasks [START_REF] Elhallaoui | Multi-phase dynamic constraint aggregation for set partitioning type problems[END_REF][START_REF] Elhallaoui | Dynamic aggregation of setpartitioning constraints in column generation[END_REF][START_REF] Benchimol | Stabilized dynamic constraint aggregation for solving set partitioning problems[END_REF], close time instants [START_REF] Macedo | Solving the vehicle routing problem with time windows and multiple routes exactly using a pseudo-polynomial model[END_REF], nearby locations, etc. However, in capacitated problems, elements with significantly different resource consumptions can hardly be considered equivalent, as the capacity constraints can substantially change the column sets covering such elements.

2IDA could simply be transformed to perform equality aggregation by imposing α j = 0, ∀j ∈ [1..k]. However, for both linearity and equality restrictions, the aggregated pricing for P k generates only constraints that are valid for P , yet P contains some solutions that cannot be projected into feasible P k solutions. Non-valid (overly-strong) P constraints implicitly arise during the P k construction by combining valid columns and exchange vectors [START_REF] Vanderbeck | A generic view of dantzig-wolfe decomposition in mixed integer programming[END_REF] that are not valid for P.

Let us illustrate the process on a valid P constraint involving some terms a i y i + a j y j , i.e., the corresponding column has the form [a i . . . a j . . . ] ⊤ (assume without generality loss that i is the first element), where the dotted elements indicates other values of the column. An equality aggregation y i = y j would makes a i and a j interchangeable in the above column. This property generates exchange vectors of the form [1 . . . -1 . . . ] ⊤ , where the dotted elements represent 0s. A valid constraint and an exchange vector produce implicit constraints as follows:

   ai . . . aj . . .    + ψ    1 . . . -1 . . .    →     ai+ψ . . . aj -ψ . . .     , (4.1)
where ψ ∈ [-a i , a j ]. The obtained constraint may be non feasible. By setting ψ to -a i or a j , one respectively replaces a i by a j , or a j by a i . Since we deal with equality restrictions, such replacements are always bidirectional (negative or positive). Using equality aggregations, the exchange vector might relax the maximum (resp. minimum) resource constraint if it replaces an element by an element of larger (resp. smaller) resource consumption.

We will now see that the dual linearity restrictions produce exchange vectors that are fractional. The constraints of (3.1) (see p. 5) can be written using the y variables only. Using Definition 1 (p. 6), we can write y j min = αw j min + β, and y j max = αw j max + β, where y j max and y j min correspond to the elements of I j of largest (respectively smallest) consumption. This gives us α = y j max -y j min w j max -w j min and β = w j min y j max -w j max y j min w j min -w j max . Then, for each element i of group j different from these two extrema, the linearity constraint is equivalent to add a constraint y i = y j max -y j min w j max -w j min w i + w j min y j max -w j max y j min w j min -w j max .

Without restricting generality, we consider that the first element of I j corresponds to y j min and the last one to y j max . Therefore, the exchange vectors produced restricted to group j have the following form:

      a j min . . . ai . . . a j max       + ψ          w j max -w i w j max -w j min . . . -1
. . .

w i -w j min w j max -w j min          →          a j min +ψ w j max -w i w j max -w j min . . . ai-ψ
. . .

a j max +ψ w i -w j min w j max -w j min          , ( 4.2)
where ψ can be positive or negative, which may respectively decrease or increase the coefficient of a i . An implicit configuration â generated this way can have fractional coefficients. Therefore, it may lead to violated patterns because it relaxes the integrity of the variables, but it cannot violate the resource constraints-i.e., C -≤ â⊤ w ≤ C + still holds. We now show this result formally. Proposition 5. Non-valid configurations â produced by linearity aggregation from a feasible configuration a have the same resource consumption as a. Such non-valid configurations cannot violate resource constraints. This is generally false for configurations produced by equality aggregations.

Proof. We first show that the equality aggregation may lead to non-valid configurations â that violate the resource constraint. It is sufficient to exhibit an example. Take C -= 8 and C + = 10, and two elements of resource consumption w 1 = 8 and w 2 = 3 in the same group. Configuration [0, 3] ⊤ is valid. A exchange vector [1, -1] ⊤ with ψ = 3 leads via (4.1) to [3, 0] ⊤ , which corresponds to a total resource consumption of 24. For the minimum resource consumption, an example can be produced by taking the valid configuration [1, 0] ⊤ and using ψ = -1 in (4.1).

We now show that non-valid configurations â generated by linearity restrictions do not violate the resource constraint because the total resource consumption cannot be modified by linearitybased exchange vectors. Recall that the configurations a generated by the pricing sub-problem are valid configurations, i.e. C -≤ a ⊤ w ≤ C + . An implicit configuration â is obtained by iteratively applying exchange vectors to a valid configuration a.

Let a = [a j min . . . a i . . . a j max ] ⊤ be a feasible configuration and e the exchange vector

w j max -wi w j max -w j min . . . , -1, . . . wi-w j min w j max -w j min ⊤
. We observe in (4.2) that combining a and e can lead to new artificial configurations the form a + ψe.

Initially the resource consumption of configuration a is a ⊤ w. Note that in a + ψe, the only coefficients to be modified are related to a j min , a j max and a i . In the original configuration, the total resource consumption of these three elements is a j min w j min + a i w i + a j max w j max . In a + ψe, this resource consumption becomes:

(a j min + w j max -w i w j max -w j min ψ)w j min + (a i -ψ)w i + (a j max + w i -w j min w j max -w j min ψ)w j max ,
which is equivalent to: a j min w j min + a i w i + a j max w j max + w j max -w i w j max -w j min ψw j min -

ψw i (w j max -w j min ) w j max -w j min + w i -w j min w j max -w j min ψw j max
This simplifies to a j min w j min + a i w i + a j max w j max , which means that the resource consumption is the same in a and a + ψe. By iteratively adding exchange vectors to the initial valid configuration a, a configuration â such that C -≤ â⊤ w = a ⊤ w ≤ C + is obtained. Therefore, any such configuration â does not violate the resource constraint.

Cutting-Stock Application and Implementation Aspects

In order to validate our methodology, we applied 2IDA to the well-known Cutting-Stock Problem. We first discuss the effectiveness of linearity restrictions compared to equality restrictions in this special case. Then, we describe in greater detail some group-splitting strategies improved by additional problem knowledge.

The cutting-stock problem is the following. We are given a bin size C ∈ N and a list 1, . . . , n of items. Each item i has a weight w i and a demand b i . The objective is to minimize the number of bins to use to accommodate all items. This problem can be solved [START_REF] Gilmore | A linear programming approach to the cutting stock problem[END_REF] by an extended formulation whose master problem is a set-covering model, and whose pricing sub-problem is a knapsack problem. Therefore, it lies in the scope of our study, fitting the initial model (2.2) very well. We use weights of the items for ressource consumption and the total capacity consumption must be between 0 and C + . To simplify the notation, we will simply note C + as C.

Linear Restrictions vs. Equality Restrictions

Since the equality aggregation is a special case of linear aggregation, the quality of the bound obtained by the latter will be better. For k = 1, the linear aggregation even has a worst-case performance that cannot be reached by the equality aggregation. Proposition 6. The asymptotic worst case performance of 2IDA with k = 1 for the standard Cutting-Stock Problem is 1/2, and this ratio is tight.

Proof. Polytope P 1 contains the solution y i = wi C . This dual solution is equivalent to the optimal solution of the model of [START_REF] Kantorovich | Mathematical methods of organizing and planning production[END_REF], which is known to have an asymptotic worst-case performance of 1/2.

To prove the tightness, consider an instance with n = 3, w = ǫ, C 2 + ǫ, C -ǫ and b = [1, M, 1], where ǫ > 0 is a sufficiently small real and M is a sufficiently large integer. The optimal solution of this instance is M + 1. We observe that the set of R 1 of feasible resource consumptions is {ǫ,

C 2 + ǫ, C 2 + 2ǫ, C -ǫ, C}. Since w 2 = C 2 +
ǫ can be as close to w1+w3 2 as desired (using a small enough ǫ), y 2 can also be arbitrarily close to y1+y3 2 -as long as we have only one group, i.e., k = 1. We now observe that y 1 + y 3 ≤ 1, based on the pattern with total consumption C. As such, we can conclude that y 2 can be arbitrarily close to 1 2 ; as such, the objective value OPT(DSCP

(P 1 )) is b ⊤ y = b 1 y 1 + b 2 y 2 + b 3 y 3 = y 1 + M y 2 + y 3 ,
which has an asymptotic value (when ǫ → 0) of M 2 + 1, while the instance has an optimum of M + 1.

This result does not hold for the equality aggregation. For k = 1, all dual values are equal to 1 ⌊C/wmin⌋ , otherwise the valid dual-constraint that uses ⌊C/w min ⌋ items of size w min would be excluded. Therefore, the ratio between the value found by the model and the optimal solution is not bounded by any constant.

Split Selection Based on Dual-Feasible Functions

In Section 3.4, we have mentioned general principles for making splitting decisions in 2IDA. For the Cutting-Stock Problem, we can also use ad-hoc methods based on the knowledge of good dual solutions (obtained by applying so-called Dual-Feasible Functions [START_REF] Lueker | Bin packing with items uniformly distributed over intervals[END_REF][START_REF] Clautiaux | A survey of dual-feasible and superadditive functions[END_REF]).

Definition 7. A function f : [0, C] → [0, 1] is a dual-feasible function (DFF), if the following holds i∈I b i w i ≤ C =⇒ i∈I b i f (w i ) ≤ 1
for any index set I such that b i ∈ Z + and w i > 0. To obtain a valid lower bound for the cutting-stock problem, it is sufficient to apply a dualfeasible function on the weights of the items and to sum up the obtained values. This way, a dual feasible function generate a valid dual solution y: it takes the weights as an input, and associate to each weight w i the corresponding dual value y i .

Given a good dual-feasible solution f , a good method for splitting groups consists of splitting [0, C] in such a way that f is linear over each the groups. Therefore, our group split method dedicated to the cutting-stock problem tries to reproduce the structure of the best dual-feasible function available. The motivation comes from the fact that most of the classical DFFs have a piece-wise linear form. The different pieces define a natural partition of the element set. By identifying all intervals of [0, C] corresponding with the different pieces of the DFF, our DFF-based split method computes a priori k intervals before starting constructing P k . If the function has k pieces, splitting the dual values into k groups leads 2IDA at least to the same feasible solutions as those returned by the DFFs. However, the best function is not always useful. For example, f (x) = x C only defines k = 1 interval. Therefore, such specific functions are not used in practice.

Split Selection Based on Outside Reference Solutions

We now introduce a split method that implements the principles from Section 3.4.2 and tunes the evaluation for Cutting-Stock. While many details below are especially tailored for Cutting-Stock, similar formulas can be used for other applications.

Since elements are sorted by increasing weight within each group, we can express a split decision by a pair (j, n j1 ), where j is the group to be split, and n j1 the index of the highest element in the first subgroup j 1 . Our splitting method assigns scores h(j, n j1 ) to all potential splits. Three types of information are used to determine a total score h(j, n j1 ) = h 1 (j, n j1 ) • h 2 (j, n j1 ) • h 3 (j, n j1 ). In the end, the potential split of maximal score is selected.

The first two scores h 1 and h 2 are based on the simple ad-hoc considerations related to the cutting-stock problem. Score h 1 (j, n ji ) is based on weight spread information, i.e., h 1 (j, n j1 ) = w j max -w j min if 1 < n j1 , n j -n j1 or 1 otherwise. The special cases of splitting a group at positions 1 or n j -1 should be avoided, since they produce a subgroup of size 1. However, in general, the probability of choosing a group j is proportional to its weight spread, as it is reasonable to split larger segments first. Score h 2 evaluates the impact of a group based on its relative position in [1..k]. The "extremal" groups j = 1 and j = k are always associated to the lightest and respectively the heaviest elements and therefore are critical in the case of cutting-stock. These two groups are given a bonus by scoring function h 2 . This part of the procedure is clearly hand-tailored for the cutting-stock and is not likely to generalize to other problems.

The major part of the score is given by h 3 . It is based on a comparison of the current optimal solution y * k of P k to an outside reference solution such that y * k → y u is an open direction. As hinted in Section 3.4.2, we assign higher scores to groups that are more likely to lead (using the process from Section 4) to implicit non-valid constraints that separate y * k from y u . For this, we compare y * k with y u . While y * k is group-wise linear structure, this is not the case for y u ; we actually try to find groups on which the difference y uy * k is far from indicating a linear transformation.

The difference between y * k and y u for a given group j is evaluated using the ∆ operator, defined by

∆ j i1,i2 (y * k , y u ) = i2 i=i1 b j i • (y * k -y u ) j i , (5.1)
where (y * ky u ) j i is the i th position of the vector y * ky u restricted to the elements of group j. To simplify the notations, we will hereafter ignore the argument (y * k , y u ) of the operator ∆ j i1,i2 ; no confusion can arise because all comparisons below concern y * k and y u . Operator ∆ is useful for detecting differences between y * k and y u and proposing useful splits. Let us analyze the case where the global ∆ evaluation over group j is positive-i.e., consider ∆ j 1,nj > 0 (the opposite case ∆ j 1,nj ≤ 0 is symmetric). This indicates a global positive "trend" on group j: the objective value of y * k could be improved by (non-linearly) increasing most dual values in group j; such operations can be performed while maintaining feasibility in P. Secondly, if ∆ j 1,i is negative for some i ∈ [1..n j ), a group split at (j, i) would allow 2IDA to decrease the first i dual values and increase the others. In this manner, 2IDA would follow the open direction y * k → y u . A more frequent situation is associated to a segment [i a , i b ] ∈ [1..n j ] such that ∆ j ia,i b < 0; to isolate dual values of [i a ..i b ] from the rest of the group, one should consider split decisions (j, i a -1) or j(j, i b ).

For this, we divide [1..n j ] in three segments k could decrease by dropping constraints that link elements of [i a ..i b ] and other elements from the group. The most fortunate case is related to i a = 1 (or, analogously, i b = n j ), because this would make [1..i a ) empty; as such, a split at i b would surely allow the first i a elements to decrease while increasing the others. The scoring function h 3 is determined by: h 3 (j, i a ) = ∆ j 1,ia-1 , h 3 (j, i b ) = ∆ j i b +1,nj and h 3 (j, i) = 0, ∀i / ∈ {i a , i b }. Finally, if the above formulas lead to i a > i b , we simply consider that the direction y * k → y u does not give any useful information, and we set

h 3 (j, i) = 1, ∀i ∈ [1..n j ].

Numerical Evaluation

We first evaluate 2IDA on the Cutting Stock Problem (CSP) in Sections 6.2-6.4. Section 6.5 continues the evaluation with experiments on the Cutting Stock Problem with Maximum Waste (CSP-MW).

We first focus on investigating the internal 2IDA dynamics (Section 6.2), i.e., we report the running time and the number of iterations (sub-problems) for the lower and the upper bounding calculations. In Section 6.3, we compare the intermediate 2IDA lower bounds with the classical Lagrangian bounds in column generation. Section 6.4 validates the fact that, in many cases, 2IDA reaches the final CG optimum (OPT CG ) more rapidly than classical Column Generation (CG).

Experimental Conditions and CSP Instances

We implemented both 2IDA and a classical non-stabilized CG algorithm based on Dynamic Programming (DP) for pricing the columns. A similar implementation approach (using lists of states) has been used for all DP routines of 2IDA, namely: the coefficient calculation preprocessing stage (Section 3.2.1), the multiple-choice aggregated pricing for optimizing DSCP(P k ) in Section 3.2.2 (or in "constraint inheritance" process from Algorithm 2), the upper bounding pricing for optimizing over P u (Section 3.5.2).

The first split decisions are made based on the intervals returned by some classical dual-feasible function (DFF, Section 5.2). We only used DFFs with relatively few intervals (maximum 15). As soon as 2IDA requires more intervals, the upper bound and the corresponding outside reference solution (when available) are used to guide further group split decisions (see Sections 3.4.2 or 5.3). We turn to classical CG if ⌈lb k ⌉ + 1 = ⌈ub k ⌉ or if k = 10. More exactly, when 10 iterations are not sufficient to reach optimality, the current set of columns used for upper bounding is used to start a CG phase that is run up to optimality.

We used a large set of instances from a dozen of benchmark sets (3127 individual instances). Their characteristics are described in great detail in Appendix A. The times we report are obtained using the same computing setting. All technical characteristics (the CPU model, LP solver, source code, etc.) are also specified in Appendix A.

2IDA Profile: Computing Effort and Usefulness of Upper Bounding

Computing an upper bound is not a mandatory part in 2IDA. Since it requires running nonaggregated column generation, it is interesting to evaluate the impact of the upper bound in term of total computing time.

Table 1 presents the computing effort in terms of CPU time and sub-problems needed, both for lower bound and intermediate upper bounding. We actually compare two 2IDA versions: one using the intermediate upper bounds from Section 3.5 ("yes" in Column 2) and a "pure" 2IDA version that does not use such intermediate upper bounds ("no" in Column 2). However, recall (Section 6.1) that we always turn to the classical CG algorithm if the optimum is not reached at k = 10. As such, the indicated upper bounding computing effort includes both the intermediate upper bounds and the upper bounds computed during the final CG process.

The first conclusion is that 2IDA is generally more efficient when intermediate upper bounds are used. These bounds are useful to guide the split decision or to stop the process more rapidly. The number of y * k -tight constraints in the polytope P u (see Section 3.5) is significantly smaller than the total number of constraints in P. Therefore, the optimization of P u does not introduce prohibitively large slowdowns. As expected, the efficiency of 2IDA is greatly improved when the optimality can be proved only using intermediate lower bounds for some k < 10.

For full convergence, 2IDA spends more computing time on upper bounding than on lower bounding. The fact that the method switches to classical CG after a number of iterations increases the difference between the UB and LB computing efforts. This difference can be significant, most notably for vb50-3 and vb50-5, where the lower bounding stops at k = 1 and classical CG is well-known specialization of the Lagrangian bound dedicated to this specific model-see Appendix B for the exact formula and references. Table 2 reports a comparison constructed by the following protocol. We first run the classical CG on each instance. This produces for each instance a reference computing time T CG and an optimal value OPT CG . Then, we run for the same instance both CG and 2IDA with a time limit of p • T CG , where p ∈ {5%, 10%, 20%, 30%, 40%}. For each such p, we note L p the best lower bound obtained after p • T CG time. If 2IDA is stopped during step k, L p is either the optimum of DSCP(P k-1 ) or the best Lagrangian dual bound obtained during step k. Table 2 reports the evolution of L p /OPT CG when p goes from 5% to 40%.

Within 5%T CG 2IDA finds the CG optimum for all instances bin3 and vb50-{2,3,4,5}. For the same computing time, the Lagrangian CG lower bounds are equal to respectively 0.89 OPT CG , 0.49 OPT CG , 0.44 OPT CG , 0.52 OPT CG , and 0.52 OPT CG on average. Within 10%T CG time, the lower bounds produced by 2IDA are better than 0.92 OPT CG for all instance sets (except m35, see below). Withing the same time limit 10%T CG , the Lagrangian GC bounds are at most 0.7 OPT CG except for instance set bin3

2IDA slowly closes the gap that is already small after 20%T CG ; this is obvious when looking at columns 20% to 40%. When the number of groups becomes larger than 10, 2IDA switches to classical column generation. The same convergence issues occur and even if the Lagrangian bound is almost equal to OPT CG , the method may need time to converge.

The only case in which the classical Lagrangian bound is better is m35, only for p = 5%. These instances have all weights larger than C 3 , and so, all columns have two non-zero coefficients. The first (almost linear) approximations of the dual values are likely to be weak. Note however that after 10% of the time needed for CG, 2IDA becomes more effective in computing dual bounds (0.80 OPT CG on average vs. 0.64 OPT CG for CG).

Computational Effort Required to Reach the Optimum

The main purpose of 2IDA is to produce high-quality dual bounds in a faster way than the classical CG. However, the fact that almost optimal dual bounds are produced in the first iterations hints that the method can help proving optimality in a faster way. converge by either method. 2IDA is significantly faster for 2402 instances out of 3127 (76%), and at least as good in 2544 cases (81%). On average, our method is faster than CG for all data sets. This can be explained by the fact that less computationally demanding (aggregated) pricing procedures are run using 2IDA. The first iterations are more effective for producing a good dual bound (as shown above). Furthermore, by only generating y * k -tight constraints during the 2IDA upper bounding process (Section 3.5), we introduce a form of implicit stabilization around those good initial solutions.

Note that the most difficult instances are hard. For these instances, our approach is always faster (a factor 3 on average) and produces less sub-problems.

To conclude, we would like to give two essential keys for the efficiency of the method. First, dual-feasible functions are useful to make good split decisions. This means that knowing a priori the structure of a good dual solution can be useful in 2IDA. Secondly, the incremental construction of polytope P k+1 from P k (Section 3.6) reduces the computing time by a wide range compared to constructing P k+1 from scratch at each iteration. This part should not be overlooked in a good implementation of 2IDA: in certain cases, this inheritance process allows 2IDA to compute lb k+1 from lb k in virtually no time (see detailed instance-by-instance results at www.lgi2a.univ-artois.fr/ ~porumbel/cs/details2ida.zip).

Cutting Stock Problem with Maximum Waste

Since 2IDA is not dedicated to the cuttig-stock problem and its structure, we now consider the Cutting Stock Problem with Maximum Waste (CSP-MW). The problem is defined as follows: all valid patterns are required to have a total length (weight) between C -and C + . In practice, one can consider that C + = C is the roll length and C + -C -is the maximum allowed waste. This constraint arises in industry when large pieces of waste cannot be recycled. The new maximum waste constraints might make certain items i be cut more than b i times (such additional cuts are not considered as waste).

It is interesting to consider this problem, because there is no equivalent of the dual-feasible functions for this variant, and thus we have no information about the structure of good dual solutions.

We consider the same instances from the previous sections. For each benchmark set, we report individual results for the first three instances. The maximum waste is set as follows: -0 for the m01 instances (i.e., C -= C + = C = 100). The instances m20 and 35 were not used because they contain only large items (at least 20 100 C or respectively 35 100 C, see Appendix A), and so, they are most often infeasible; -2 for all vb50 instances (i.e., C -= 99998); -0.5%C for vbInd instances; for this set, we provide results for the only three instances with n ≥ 30; -4, for the hard instances (i.e., C -= 196).

Table 4 compares the bounds obtained by 2IDA without upper bounding and with splitting decisions taken using weight spread indicators (as described in Section 3.4.1). The conclusions are clear: the 2IDA bounds clearly dominate the Lagrangian bounds. More exactly, even the first 2IDA bounds for k = 1 (Column 2) are usually larger than the Lagrangian bounds reported after 2-3 times more time (Column 7). The only exception arises for the hard instances, where CG converges more rapidly.

Conclusions and Perspectives

We described an aggregation method for computing dual bounds in column generation in a way that is faster than classical Lagrangian bounds. The new method proceeds by constructing an iterative inner approximation of the dual polytope and converges towards the optimum of the original model. Several computational techniques are used to incrementally construct these inner polytopes and obtain a practically effective method. Computational results show that, besides generating fast dual bounds, the method can reach the optimum of the original model more rapidly than classical column generation in many cases.

Future research work will concentrate on the case with several resources, and also more general cases where the pricing sub-problems are solved by an ILP solver.

VB-a instances vb10, vb20 (random instances [START_REF] Vanderbeck | Computational study of a column generation algorithm for bin packing and cutting stock problems[END_REF]) and vbIndst (industrial instances [START_REF] Vanderbeck | Computational study of a column generation algorithm for bin packing and cutting stock problems[END_REF]). We selected only the industrial instance with one bin type;

VB-b instances vb50-1, vb50-2,. . . , vb50-5, random instances;

BP three randomly generated bin-packing instance sets [START_REF] Clautiaux | A survey of dual-feasible and superadditive functions[END_REF]: m01, m20, m35; hard 10 bin-packing instances [START_REF] Scholl | BISON: a fast hybrid procedure for exactly solving the one-dimensional bin-packing problem[END_REF] where: λ is the primal variable vector of the column generation, and M RC = min a∈A 1a ⊤ y ≤ 0 is the non-positive minimum reduced cost at the current iteration, The key for calculating the Farley bound resides in observing that the inequality below can be artificially added to the initial LP.

  a ≥ b i , ∀i ∈ 1, . . . , n λ a ∈ N,

Algorithm 2 :

 2 Two-Phase CG for optimizing DSCP(P k+1 ): construct P ′ k+1 then P k+1 Data: Optimal DSCP(P k ) solution [α 1 . . . α k β 1 . . . β k ] ⊤ , existing constraints A base k Result: OPT(DSCP(P k+1 )) Phase 1: Lift current solution and describe P ′ k+1 Express current [α 1 . . . α k β 1 . . . β k ] ⊤ as a feasible DSCP(P k+1 ) solution:

  [1..i a ), [i a ..i b ], (i b ..n j ], determining i a and i b as follows:i a = max{i ∈ [1..n j + 1] : ∆ j 1,i-1 ≥ 0} and i b = min{i ∈ [0..n j ] : ∆ j i+1,nj ≥ 0}. Let us first consider i a ≤ i b . The middle segment [i a ,i b ] has an "inverse trend" compared to [1..n j ], i.e., y * k could feasibly improve in P, by decreasing its values over [i a , i b ], while increasing its values over the rest of [1..n j ]. Roughly speaking, we are looking for a sequence of indices [i a , ..i b ] on which y *

λ≥0 a∈A ( 1

 1 a ⊤ y)λ a + b ⊤ y ≥ min λ≥0 M RC a∈A λ a + b ⊤ y, (B.1)

  a∈A λ a ≤ OPT CG , (B.2) The 2IDA Lagrangian bounds from Section 6.3 are calculated in the same manner, i.e., the whole reasoning still holds by only using A k instead of A. Combining (B.2) with M RC ≤ 0, we can reduce (B.1) to OPT CG ≥ M RC • OPT CG + b ⊤ y. The Farley lower bound L F y follows immediately: L F y = b ⊤ y 1-MRC ≤ OPT CG .

Table 2 :

 2 Comparison of intermediate dual bounds for CG an 2IDA. We report the average ratio L p /OPT CG for each set of instances for p =5% to p =50% of the time T CG . For example, row vb20 of column "5%" shows that if both methods are run on instances of data set vb20, and stopped after 5% • T CG , CG will output a dual bound equal to 0.41 OPT CG , while 2IDA will output a dual bound equal to 0.64 OPT CG .

			Percentage of CG time T CG needed for full convergence
		5%	10%	20%	30%	40%
	Instance set CG 2IDA CG 2IDA CG 2IDA CG 2IDA CG 2IDA
	vb20	0.41 0.64 0.55 0.95 0.70 0.99 0.85 0.99 0.94 0.99
	vb50-1	0.29 0.88 0.42 0.97 0.65 0.98 0.79 0.98 0.84 0.99
	vb50-2	0.49 1.00 0.70 1.00 0.87 1.00 0.97 1.00 0.98 1.00
	vb50-3	0.44 1.00 0.57 1.00 0.65 1.00 0.88 1.00 0.98 1.00
	vb50-4	0.52 1.00 0.66 1.00 0.76 1.00 0.94 1.00 0.98 1.00
	vb50-5	0.52 1.00 0.60 1.00 0.65 1.00 0.79 1.00 0.97 1.00
	vbInd	0.52 0.64 0.70 0.87 0.71 0.93 0.82 0.93 0.88 0.99
	m01	0.41 0.94 0.51 0.95 0.65 0.97 0.73 0.97 0.77 0.98
	m20	0.49 0.74 0.55 0.92 0.66 0.95 0.75 0.98 0.79 0.98
	m35	0.64 0.47 0.64 0.80 0.69 0.92 0.73 0.96 0.75 0.98
	hard	0.89 1.00 0.89 1.00 0.89 1.00 0.89 1.00 0.94 1.00

Table 3

 3

		Avg. CPU [ms] 2IDA vs CG (CPU) nb sub-problems
	Instance set 2IDA	CG	≺	≃	≻	2IDA	CG
	vb20	9161	11653 10	4	11	89	63
	vbIndstr	7580	14252 10	1	6	43	51
	vb50-1	39827	50164 14	2	4	179	125
	vb50-2	89432	125399 17	0	3	206	223
	vb50-3	44068	106356 20	0	0	80	159
	vb50-4	66498	87604 19	1	0	242	199
	vb50-5	28490	85680 20	0	0	62	160
	m01	818	849 549 43	408	161	191
	m20	423	620 793 49	158	128	173
	m35	207	379 920 42	38	71	134
	hard	117908 376112 10	0	0	219	768

compares the CPU time measures and the number of sub-problems calls needed to

Table 3 :

 3 The computing effort for final convergence with 2IDA and CG. Columns "Avg. CPU" report the time (in ms.). Columns "2IDA vs CG (CPU)" report the number of times 2IDA is (1) faster than CG (column ≺), (2) roughly equally fast (difference < 5%, column ≃), and (3) slower (column ≻). The last two columns compares the total number of pricing calls for 2IDA and CG. For 2IDA, the total number of pricing calls considers both Multiple-Choice Knapsack problems (columns of P k ) and classical Knapsack problems (columns used for upper bounding, including the inherited columns from Section 3.6).

Table 4 :

 4 The first column indicates the instance name. Columns 2, 4 and 6 respectively provide the bound obtained by 2IDA for k from 1 to 3. Columns 3, 4, and 5 respectively report the Lagrangian bound obtained using a similar CPU time as 2IDA. More exactly, for each k, we report the best Lagrangian bound found during CG at the first moment when the CG running time is larger than the one needed by 2IDA. The last two columns indicate the CG optimum with and respectively without maximum waste. * Remark that these 2IDA bounds are better than the CG optimum in the last column. Any bound for the standard Cutting Stock Problem (e.g., based on DFF) could never reach this value, as they would all be limited by the Cutting-Stock optimum.

Table 5 :

 5 long-acknowledged for their difficulty in the bin-packing literature General characteristics of the instances. Columns 4 and 5 indicate the (approximate) interval of the demand value, and, respectively, item weights. In the archive http://www.math.u-bordeaux1.fr/~fvanderb/data/randomCSPinstances.tar.ZB The Lagrangian Bound for CSP and CSP-MWWe now recall the specialized Lagrangian lower bound used for the standard Cutting-Stock Problem and that can be used for the Cutting-Stock Problem with Maximum Waste. We come back to the main Set-Covering column generation model (2.1)-(2.2), p. 3. Similarly to what is done in [2, § 2.2], [26, §. 3.2], [17, §. 2.1] or [4, § 1.2], Lagrangian bounds are constructed as follows. Consider a given iteration of the CG process with dual values y. We use these values y as multipliers of the Lagrangian relaxation of (2.1). The relation between the CG optimum OPT CG and the Lagrangian bound L y can be written: OPT CG ≥ L y = min

	Name	n	C	avg. b span avg. w span	Description
	vb10	10 10000	[10, 100]	[1, 1 2 C]	20 random instances [26]: CSTR10b50c[1-5]* files a
	vb20	20 10000	[10, 100]	[1, 1 2 C]	25 random instances [26]: CSTR20b50c[1-5]* files a
	vb50-c1	50 10000	[50, 100]	[1, 3 4 C]	20 random instances [26]: CSTR50b50c1* files a
	vb50-c2	50 10000	[50, 100]	[1, 1 2 C]	20 random instances [26]: CSTR50b50c2* files a
	vb50-c3	50 10000	[50, 100]	[1, 1 4 C]	20 random instances [26]: CSTR50b50c3* files a
	vb50-c4	50 10000	[50, 100]	[ 1 10 C, 1 2 C]	20 random instances [26]: CSTR50b50c4* files a
	vb50-c5	50 10000	[50, 100]	[ 1 10 C, 1 4 C]	20 random instances [26]: CSTR50b50c5* files a
	vb50-b100 50 10000	[1, 210]	[ 1 10 C, 1 2 C]	20 random instances [26]: CSTR50b100c4* files a
	m01	100 100	1	[1, C]	1000 random bin-packing instances [7];
	m20	100 100	1	[ 20 100 C, C]	1000 random bin-packing instances [7];
	m35	100 100	1	[ 35 100 C, C]	1000 random bin-packing instances [7];
	Hard	≈ 200 100000	1 -3	[ 20 100 C, 35 100 C]	Bin-packing instances, known to be more difficult

a

The second column indicates whether the upper bound from Section 3.5 is used to guide the search, i.e., for splitting and for early stopping. Columns 3-5 provide the total CPU time used for lower bounding and upper bounding calculations. The last 6 column present the minimum, average and maximum number of pricing calls in each part of 2IDA. Recall that UB pricing calls corresponds to both Section 3.5 and the classical CG steps that are run at the end of the method.

run immediately because ⌈lb k ⌉ + 1 = ⌈ub k ⌉. The running time dedicated to upper bounding is almost ten times larger than the lower bounding time. The construction of P 1 is not carried out by column generation, since the model has only 2 variables and it is faster to generate directly all constraints. This explains why the number of calls to multichoice knapsack routines is sometimes 0 (columns "LB (2IDA) calls"). In such cases, most of the total 2IDA computing effort is actually spent on proving that this lower bound is optimal. The 2IDA version with no intermediate upper bounds may be efficient when small instances are considered. The smallest instances vb10 with n = 10 could be solved by using no upper bounding at all. In this case, 2IDA reaches rapidly a state in which all groups have at maximum two elements and the 2IDA lower bound is optimal. This also happens for the smaller vbInd instances. Clearly, the results on the most difficult instances are more useful. Therefore, in the remaining we only focus on the version with intermediate upper bounds and we ignore the small instances vb10.

Quality of the Intermediate Dual Bounds

We compare the dual bounds iteratively produced by 2IDA with the Lagrangian dual bounds produced during CG each time the sub-problem is solved. We used the Farley bound, which is a

A Instance Sets and Computing Setting

All reported CPU times are obtained using the following computing environment: a HP ProBook 4720 laptop clocked at 2.27GHz (Intel Core i3), with the gnu g++ C++ compiler on Linux Ubuntu, kernel version 2.6. The master problems are solved using Cplex 12.3 and Ilog Concert libraries. The source code and the instances are publicly available on-line at www.lgi2a.univ-artois.fr/ ~porumbel/cs/. Detailed results (instance by instance) are also publicly available at this address.

We now describe in details the instances we used in this paper. They form a set of 3127 individual instances. The number of items ranges from 10 to 200 and the capacities is between 100 and 100000.