Column generation decomposition with the degenerate constraints in the subproblem

European Journal of Operational Research 207 (2010) 37–44

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Continuous Optimization

Column generation decomposition with the degenerate constraints in the subproblem Abdelmoutalib Metrane a, François Soumis b, Issmail Elhallaoui b,* a b

Université Hassan 1, Ecole Nationale des Sciences Appliquées de Khouribga, Hay El Jadid, El Biout ‘‘OCP”, B.P. 77, Khouribga, Morocco École Polytechnique de Montréal and GERAD, Département de Mathématiques et de Génie Industriel, C.P. 6079, Succ. Centre-Ville, Montréal, Québec, Canada H3C 3A7

a r t i c l e

i n f o

Article history: Received 27 April 2009 Accepted 1 May 2010 Available online 13 May 2010 Keywords: Linear programming Primal simplex algorithm Column generation Degeneracy

a b s t r a c t In this paper, we propose a new Dantzig–Wolfe decomposition for degenerate linear programs with the non degenerate constraints in the master problem and the degenerate ones in the subproblem. We propose three algorithms. The first one, where some set of variables of the original problem are added to the master problem, corresponds to the Improved Primal Simplex algorithm (IPS) presented recently by Elhallaoui et al. [7]. In the second one, some extreme points of the subproblem are added as columns in the master problem. The third algorithm is a mixed implementation that adds some original variables and some extreme points of a subproblem to the master problem. Experimental results on some degenerate instances show that the proposed algorithms yield computational times that are reduced by an average factor ranging from 3.32 to 13.16 compared to the primal simplex of CPLEX. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction The primal simplex algorithm introduced by Dantzig in 1947 (see [2]) is often used for solving linear programs. It is, in particular, the most efficient algorithm for re-optimizing after adding cuts in a Branch and Bound schema or new columns in a column generation method. This algorithm iteratively improves the objective value by moving from basis to basis. The time per iteration grows with the number of variables and even more with the number of constraints. However, the primal simplex may execute many iterations without improving the objective value when the basic solution is degenerate. To improve the solution time, research has focused first on reducing the number of variables to take into account at each iteration. Partial pricing is used to reduce the pricing time. Column generation [8,3] uses a master problem restricted to a subset of variables and one or many subproblems to identify variables to add to the master problem. Research on reducing the number of constraints to take into account at each iteration is more recent even though this factor has more impact on the solution time. Three papers by Elhallaoui et al. [4–6] recently propose many dynamic constraint aggregation methods reducing greatly the solution time of set partitioning problems. Most recently Elhallaoui et al. [7] presented the

* Corresponding author. E-mail addresses: [email protected] (A. Metrane), [email protected] (F. Soumis), [email protected], [email protected] (I. Elhallaoui). 0377-2217/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2010.05.002

Improved Primal Simplex (IPS) for general linear programming problems combining ideas from Elhallaoui et al. [5] and Pan [12]. The authors in [7] claimed that their straightforward implementation of IPS has allowed reduction of the solution time by a factor of 3 compared to the primal simplex of CPLEX. Raymond et al. [13] proposed an improved version of IPS allowing higher reduction factors. It is this improved version that will be compared to the algorithms we are proposing in this paper. The idea behind IPS is to remove all or a part of the degenerate variables from the basis. Let p 6 m be the number of remaining basic variables. A reduced basis is created with the p constraints that remain linearly independent after removing these degenerate variables. A variable is said to be compatible with a reduced basis if it is linearly dependent on the columns of this reduced basis. In other words, it can be pivoted into the reduced basis and take a positive value without violating any of the removed constraints. Otherwise, it is said to be incompatible. In IPS, we have two kinds of iterations. It executes minor iterations when solving the reduced problem containing only the p constraints and the compatible variables by using the simplex algorithm. It performs major iterations in order to identify a subset of incompatible variables allowing an improvement of the objective value of the reduced problem. To do so, a complementary problem, containing m  p + 1 constraints and all incompatible variables, is solved. The identified incompatible variables and also some necessary constraints making these variables compatible are added to the reduced problem. A solution is optimal for the original problem when no incompatible variables improving the objective value can be identified by the complementary problem. Note also

38

A. Metrane et al. / European Journal of Operational Research 207 (2010) 37–44

that during a major iteration, some degenerate variables and constraints could be switched from the reduced problem to the complementary problem. Numerical experiments show that the solution time of the vehicle and crew scheduling and fleet assignment instances tested is reduced by impressive factors ranging from 4 to 10 even though this general method does not exploit the underlying structure of the test problems. In this paper, we show that the IPS algorithm can be seen as a special case of column generation methods in the sense of Löbel [11] and De Carvalho [14]. In fact, we can see the reduced problem as the master problem and the complementary problem as the subproblem. We also show that these authors use a master problem containing the variables of the original problem rather than the extreme points of the subproblems. For vehicle routing problems, Löbel [11] uses a master problem where the variables are arcs and the subproblem generates paths with negative reduced costs. They then add the arcs of these paths to the master problem. The authors observe that the solution process is faster when using the original variables in the master problem than using the extreme points generated by the subproblems. This formulation of the master problem necessitates that all constraints of the original problem be present in the master problem. This technique is efficient only when the subproblem has few constrains (as in knapsack or shortest path problem) and the resulting master problem does not have too many constraints. In this paper, we show that IPS is a Dantzig–Wolfe decomposition where the original problem is partitioned, in a natural way, according to its inherent algebraic structure rather than by the modeler. Furthermore, this decomposition is updated dynamically during the solution process because the algebraic structure changes as the basis changes. In this decomposition, we form a less or non degenerate master problem, containing only non degenerate constraints, that can be solved by the primal simplex algorithm and a degenerate subproblem that can be solved by the dual simplex or other specialized algorithms. We present also a variant of IPS where the extreme points of the subproblem are added as variables to the master problem. This classical implementation of Dantzig–Wolfe decomposition is more efficient for a certain class of instances than the original IPS algorithm. Finally, this paper describes a hybrid implementation that adds to the master problem some mixture of original variables and extreme points of the subproblem. We use a simple criterion to decide what kind of variables (original or extreme point) to add to the master problem. A brief description of this criterion is given in Section 6.2. This hybrid version induces a further improvement of the solution time compared to the versions using only one kind of variables. This paper permits using all the expertise developed in the last decades to reduce the number of variables for large scale linear programs in order to use it again to reduce the number of constraints in degenerate linear problems.

 x ¼ B1 b. Without loss of generality, assume that there are p positive-valued basic variables in this solution indexed from 1 to p. This implies that there are m  p degenerate basic variables (indexed from p+ 1 to m) and, as usual, n  m nonbasic variables. Q PB Let B1 ¼ where the submatrix Q PB is composed of the Q ZB 1 first p rows of B , that is, those corresponding to the positive-valued basic variables ( xj > 0; j ¼ 1; . . . ; p). This submatrix is used to define the constraints of the reduced problem, while the submatrix Q ZB is used to determine the variables involved in this problem, namely, those that are compatible with Q ZB according to the following definition. Definition 2.1. Let M be an m0  m matrix such that m0 < m. A column Aj of A is said to be compatible with M if MAj = 0. The corresponding variable xj is also said to be compatible with M. M is called a compatibility matrix. It is easy to see that the positive-valued basic variables are compatible with the matrix M ¼ Q ZB while the degenerate ones are incompatible. One can also observe that a column Aj of A is compatible with Q ZB if and only if it is linearly dependent on its first p columns A‘, ‘ = 1, . . . , p. To distinguish between the variables compatible and incompatible with Q ZB , we adopt the following notation. The vector of compatible (resp. incompatible) variables is denoted by xC B (resp. xIB ) and its cost vector by cC B (resp. cIB ). The coefficient matrix A is subdivided into two submatrices AC B and AIB such that the columns of AC B are compatible with Q ZB while the columns of AIB are incompatible with Q ZB . To obtain the reduced problem for the basis B, we first transform the constraint set (2.2) as follows:

ð AC B

AI B Þ



xC B xI B

min

xC ;xI B

c>CB xC B þ c>IB xIB ;

B

s:t:

Q PB AC B xCB þ Q PB AIB xIB ¼ Q PB b; Q ZB AIB xIB ¼ 0; xCB ; xIB P 0:

The reduced problem, denoted by RPB, is deduced from this last formulation by simply omitting the incompatible variables and removing the second set of constraints: > ðRP B Þ zRP B ¼ min cC B xC B ;

xC B P 0:

Let (LP) be a linear program defined as follows:

c> x;

s:t:

ð2:1Þ

Ax ¼ b;

ð2:2Þ

x P 0; n

n

m

ð2:3Þ mn

where x 2 R ; c 2 R ; b 2 R ; A 2 R (aij is an element of A) and x is the vector of variables. The dual variable vector associated with the constraint set (2.2) is denoted by p. The IPS algorithm uses a reduced problem which coincides with (LP) when there is no degeneracy. This problem is defined according to a feasible basis B that provides the basic feasible solution

ð2:4Þ

B

s:t: Q PB ACB xC B ¼ Q PB b;

2. The reduced problem

x

    Q PB x ¼ b () ð AC B AIB Þ C B Q ZB xI B    Q PB Q P B A C B xC B þ Q P B A I B xI B ¼ Q P B b ¼ b () Q ZB Q ZB AC B xC B þ Q ZB AIB xIB ¼ Q ZB b:

Since the columns of AC B are compatible with Q ZB , we find that Q ZB AC B ¼ 0. Furthermore, since  xj ¼ 0, for j = p + 1, . . . , m, we get that Q ZB b ¼ 0. Hence, (LP) is equivalent to:

xC

ðLPÞ zLP ¼ min



ð2:5Þ ð2:6Þ

Thus, RPB only involves compatible variables. These variables have the property that they can be entered into a basis of RPB and take a positive value without violating any of the removed constraints. This is not the case for the incompatible variables. 3. The complementary problem Let xC B be an optimal basic feasible solution for the reduced problem RPB, B* the corresponding optimal basis, and p* the number of positive-valued variables in this solution (p* 6 p). Assume that these p* variables are the first p* basic variables associated

A. Metrane et al. / European Journal of Operational Research 207 (2010) 37–44

with the columns A‘, ‘ = 1, . . . , p*. As in the preceding section, the   Q PB , where matrix (B*)1 is partitioned into two parts: ðB Þ1 ¼ Q ZB the submatrix Q PB is composed of the first p* rows of (B*)1. The columns Aj of A that are linearly dependent on the first p* columns of A are the columns compatible with the p*  m matrix Q ZB . Denote by P B and Z B the index sets of those variables that take a positive value and a zero value in xC B , respectively. Denote also by I B the index set of the variables incompatible with Q ZB . Then xC B completed by zero is also optimal for (LP) if and only if all reduced costs, i.e., corresponding to all variables, compatible or not, are nonnegative. In other words, there must exist dual variables p such that

cj  ATj p ¼ 0 for all j 2 P B ;

ð3:1aÞ

cj  ATj p P 0 for all j 2 Z B ;

ð3:1bÞ

cj 

ATj

p P 0 for all j 2 I B :

ð3:1cÞ



From [7, Proposition 3.1], if p satisfies (3.1a), it also satisfies (3.1b). The second set of conditions thus becomes redundant and may be removed. The authors of [7] propose the following complementary problem, denoted by CP B :

ðCPB Þ zCP B ¼ max y;

ð3:2Þ

y;p

>

s:t: cj  p Aj ¼ 0; cj  p> Aj P y;

8j 2 P B ; 8j 2 I B  :

ð3:3Þ ð3:4Þ

where p is the dual variables vector, cj is the cost of xj and Aj its associated column in matrix A. The following problem is the simplified dual of CP B denoted SDB :

  1  ðSDB Þ minv c>I  c>P AP1 AI1 v ;     1 AI1  AI2 v ¼ 0; s:t: AP2 AP1 e> v ¼ 1;

v 2 RjI

B j

is

the

ð3:6Þ ð3:7Þ

v P 0; : where

ð3:5Þ

ð3:8Þ dual

variable of constraint (3.4),  P A1 e> = (1, . . . , 1), cP ¼ ðcj Þj2P B ; cI ¼ ðcj Þj2I B ; ¼ ðaij Þ16i6m where AP2 j2P B  I A1 P ¼ ðaij Þ16i6m . A1 is a p*  p* invertible matrix, and similarly AI2 j2I B Observation: Let v* be an optimal basic solution to SDB . This solution is a convex combination of incompatible columns which is linearly dependent on the positive-valued columns. This convex combination AI v is thus compatible with Q ZB [7]. In fact, it is a convex combination according to the constraint (3.7). The constraints (3.6) ensure that this convex combination has coefficients with zero value in the degenerate constraints. In other words, this convex combination is in the subspace generated by the basis of the reduced problem and hence compatible with Q ZB . As the optimal values of CP B and SDB are equal, it maximizes the minimum reduced cost of the incompatible variables w.r.t. all bases of (LP) containing the positive-valued variables xj ; j 2 P B . The following notation is adopted throughout this paper:

XB ¼ fx 2 Xjx is compatible with Q ZB and cx < 0g; n o P P P withX ¼ x ¼ j2I B v j Aj jv j > 0 and j2I B v j ¼ 1 ; cx ¼ j2I B v j cj and cj is the reduced cost of incompatible column Aj. An element of XB (if exists) is a convex combination of incompatible columns compatible with Q ZB and has a negative reduced cost. It is easy to see that XB – £ if and only if zCP B < 0.

39

Algorithm 1: Improved Primal Simplex 1: Find an initial basis B. 2: Construct and solve the reduced problem RPB and update XB . 3: If XB ¼ £, then stop: the current solution is an optimal solution to (LP). 4: Choose one x of XB and build a new basis B0 from the incompatible columns Aj that form x. 5: B: = B0 . Go to Step 2

Theorem 3.1. [[7]] (a) Assume that zCP B < 0. Let v* be an optimal basic feasible solution to SDB and S B be the subset of the indices j 2 I B for which v j > 0. The set of columns Aj ; j 2 S B , is minimal in the sense that no convex combination of a strict subset of its columns is compatible with Q ZB . (b) Let x* be a non-optimal basic feasible solution to (LP) (with corresponding basis B*) obtained at an iteration of the primal simplex algorithm. Assume that x* has p* 6 m positive-valued variables, where m is the number of constraints in (LP). Entering at most m  p* + 1 variables into the current basis B* is sufficient to decrease the objective value or prove that (LP) is unbounded. Furthermore, these variables must be such that: (i) a convex combination of their corresponding columns in A must be linearly dependent on the columns of the p* positivevalued variables and ii) the maximum reduced cost of the convex combinations of variables w.r.t. all bases containing the p* positive-valued variables is negative. 4. Relationship between IPS and column generation decomposition In this paper, we show that IPS is equivalent to a column generation decomposition where the original problem is naturally partitioned according to its inherent algebraic structure. This decomposition is dynamic in the sense that it is updated automatically when the algebraic structure of the problem changes. In this decomposition, we form a less or non degenerate master problem called here a reduced problem that contains non degenerate constraints and can be solved by the primal simplex algorithm; and a degenerate subproblem (the complementary problem) that can be solved by the dual simplex or other specialized algorithms. At each iteration, some of the original variables (incompatible) are added to the master problem in a non-aggregate way. We also present a variant of IPS where the extreme points of the subproblem are added as aggregate columns to the master problem. This aggregated implementation of IPS corresponding to the Dantzig– Wolfe decomposition is more efficient for a certain class of instances than the original IPS algorithm. 4.1. Column generation decomposition with non-aggregate columns In the literature, column generation and Dantzig–Wolfe decomposition are generally decided by the modeler according to sub structure appearing in the model. In contrast to that, the IPS algorithm can be seen as a decomposition where the problem (LP) is decomposed according to its inherent algebraic structure. Actually, in Dantzig–Wolfe decomposition, the problem is decomposed into a master problem and a subproblem. By applying a Dantzig–Wolfe decomposition to a linear transformation of (LP), we show in this subsection that the resulting subproblem is equivalent to the complementary problem presented in the previous section.

40

A. Metrane et al. / European Journal of Operational Research 207 (2010) 37–44

In fact, the problem (LP) can be rewritten as follows:

min c>P xP x ;x P

I

AP1 xP AP2 xP

s:t:

þ

c>Z xZ

þ AZ1 xZ þ AZ2 xZ

þ

c>I xI ;

þ AI1 xI ¼ bP ; þ AI2 xI ¼ bN ;

xP ; xZ ; xI P 0; where

cZ ¼ ðcj Þj2ZB ; xP ¼ ðxj Þj2P B ; xZ ¼ ðxj Þj2Z B ; xI ¼ ðxj Þj2I B ; bP  Z A1 ¼ ðbi Þi2f1;::;pg ; bN ¼ ðbi Þi2fpþ1;::;mg ; ¼ ðaij Þ 1 6 i 6 m . AZ2 j 2 Z B Provided that AP1 is invertible, we can state that:

min c>P xP þ c>Z xZ þ c>I xI ; xP ;xI

s:t: AP1 xP þ AZ1 xZ þ AI1 xI ¼ bP ;     AP2 ðAP1 Þ1 AZ1  AZ2 xZ þ AP2 ðAP1 Þ1 AI1  AI2 xI ¼ AP2 ðAP1 Þ1 bP  bN ; xP ; xZ ; xI P 0: Since ðxP ; 0Þ is solution to (LP), we have: AP1 xP ¼ bP ; AP2 xP ¼ bN . This implies that bN  AP2 ðAP1 Þ1 bP ¼ 0.   Remark 4.1. It is interesting to see that AP2 ðAP1 Þ1 ; Imp is a     AP1 compatibility matrix. In fact, AP2 ðAP1 Þ1 ; Imp ¼ 0. So, AP2 every positive-valued column is compatible, and then all columns that depend linearly on these positive-valued columns are also compatible.  According to Remark 4.1, we can conclude that AP2 ðAP1 Þ1 AZ1  AZ2 ¼ 0. Thus, (LP) is finally equivalent to:

c>C xC þ c>I xI ;

min xC ;xI

s:t:

AC1 xC þ AI1 xI ¼ bP ;     1 AI1  AI2 xI ¼ 0; AP2 AP1 xC ; xI P 0;

ð4:1Þ ð4:2Þ ð4:3Þ ð4:4Þ

where the index C refers to compatible columns, i.e., columns with the index j in P [ Z. One can solve this problem using Dantzig–Wolfe decomposition. For instance, we can put the constraints set (4.3) in the subproblem and of course the first set (4.2) will be considered in the master problem. This way the subproblem becomes:

min ðc>C  k> AC1 ÞxC þ ðc>I  k> AI1 ÞxI ; xC ;xI

s:t:

ðAP2 ðAP1 Þ1 AI1  AI2 ÞxI ¼ 0; xC ; xI P 0;

where k is the dual vector associated with the master problem constraints. It is easy to see that the compatible columns have nonnegative reduced costs when we solve to optimality the master problem. This means that we can simply remove these compatible columns from the subproblem. Hence, the subproblem becomes:

min xI

s:t:

ðc>I  k> AI1 ÞxI ;   AP2 ðAP1 Þ1 AI1  AI2 xI ¼ 0; xI P 0: AP1

P 1 As is a basis of the reduced problem, then k> ¼ c> is P ðA1 Þ the dual vector associated with this basis. By substitution, we obtain the objective (3.5) and the constraints (3.6) and (3.8) of the

complementary problem presented in the previous section. So, the only difference between this subproblem and the complementary problem is the presence of the convex combination constraint (3.7) in the latter one. Observe that the subproblem is a cone. So, this convex combination constraint can be seen as a normalization that transforms each extreme ray of the subproblem to an extreme point. By adding the convex combination constraint to the subproblem, the latter becomes bounded. This concludes that this subproblem and the complementary problem are equivalent. By abuse of language, we refer hereunder to an extreme ray of the subproblem as an extreme point. In this decomposition, the master problem is the reduced problem RPB with less or no degeneracy on which the primal simplex is efficient, and the subproblem is the complementary problem SDB . We should emphasize here that the master problem is formulated according to the original variables as in [11,14]. The subproblem SDB generates an extreme point by minimizing the reduced cost. This extreme point is not added to the master problem but it is the original incompatible columns forming this extreme point that are added instead. As the subproblem is a cone, there is no need to add the convex combination constraints to the master problem, that is indeed the reduced problem. The subproblem is degenerate; it can be solved by the dual simplex or other specialized algorithms. By solving the subproblem, we either identify a set of incompatible columns that, when added to the master problem, decrease the objective function or prove optimality; in this case, the optimal solution of the master problem, completed by zeros, is an optimal solution to (LP). Note also that some constraints, making these variables compatible, are added to the master problem. The IPS will be referred to as IPS-N, which stands for IPS-Non-aggregated. The master problem and the subproblem are adjusted dynamically at each iteration. Actually, constraints and variables move from the master problem to the subproblem and vice versa during the solution process. Practically speaking, updating the decomposition at an iteration is costly. So as a trade-off, we keep some degenerate variables (or constraints) in the master problem as long as their presence is not harming the solution time due to degeneracy. However, the extreme point generated by the subproblem can be entered into the master problem as an aggregate column. In this case, we do not need to add constraints to the master problem because the aggregate column is compatible. This issue is the main topic of the next subsection. 4.2. Column generation decomposition with aggregate columns In this subsection we present a new variant of IPS that can be seen as classical Dantzig–Wolfe decomposition where the columns that we add to the master problem RPB are the extreme points of the subproblem SBB . We will refer to this variant as IPS-A, which stands for IPS-Aggregated. The algorithm pseudo-code is given in Algorithm 2. The infeasible and unbounded cases have been left out to simplify the presentation of this algorithm. They can, however, be detected in Step 1 for the infeasible case and in Steps 2 for the unbounded case. Steps 1–5 of the IPS-A algorithm are briefly discussed below. Algorithm 2: Aggregated Improved Primal Simplex 1: Find an initial feasible basis B and construct RPB and SDB . 2: Solve the reduced problem RPB. 3: Update the reduced costs of the incompatible columns and update XB . 4: If XB ¼ £, then stop: the current solution is an optimal solution to (LP). 5: Add one or many columns x of XB to RPB. Go to Step 2.

41

A. Metrane et al. / European Journal of Operational Research 207 (2010) 37–44

Step 1: The algorithm starts by finding a feasible basis B for (LP). This basis can be obtained by performing a phase I with the primal simplex algorithm. Once B is obtained, the reduced problem RPB is then constructed by defining the constraint set (2.5) and, identifying the incompatible variables and removing them. The complementary problem SDB is built as described in Section 4. Step 2: In this step, RPB is solved using the primal simplex algorithm to yield an optimal basis B* and an optimal solution xC B . Step 3: When we solve the reduced problem RPB, we obtain new dual variables. The reduced costs of the incompatible columns are updated; then only the objective function of SDB is changed. Consequently, the set XB is updated accordingly. Observe that the constraints (3.6)–(3.8) are unchanged. So, B* is identical to B opposite to the Algorithm 1. Steps 4–5: If zSDB P 0 then XB ¼ ; and the current solution is an optimal solution to (LP). Otherwise, we construct the P aggregate column xv  ¼ j2SB v j Aj from the optimal basic solution v* of SDB . We add this new column to RPB. We do not add any constraints to RPB since this column is compatible. The complementary problem SDB is solved using the dual simplex algorithm. 4.3. Accelerating strategy Preliminary tests have shown that, as in a column generation method, many iterations of the IPS-A and IPS-N algorithms may be needed to reach optimality for certain instances, yielding high computational times. To reduce the number of iterations and speed-up IPS-A and IPS-N, we propose an adaptation of a well known column generation acceleration strategy that simply consists of generating more than one column x of XB per iteration. To obtain many columns, the complementary problem is solved recursively after adding the additional constraints forbidding previous solution. More precisely, let C be the cumulative index set of the nonzero components of the solutions to the modified SDB and let d > 0 be a parameter. As long as jCj 6 d, we solve a modified SDB by adding the constraints vi 6 0 for all i 2 C to SDB . If feasible and its objective function value is negative, we construct an aggregate column xv as in Step 4. At the end of this procedure, we get one or more aggregate columns in Step 4 of Algorithm 2. These columns will be added to the reduced problem. In IPS-N, we apply the same acceleration strategy described above. The only difference is that in IPS-N, we add the incompatible columns forming these aggregate columns and build a new basis B0 in Step 4.

new problem RP AB , we obtain a new objective function value denoted by zARP . The following lemma shows that we cannot find another convex combination of the incompatible columns Aj ; j 2 S B which is compatible with Q ZB . We show then that the objective values at the first iteration of IPS-N and IPS-A are equal. Lemma 5.1. The only compatible column formed by a convex combination of incompatible columns Aj ; j 2 S B is xv  : P Proof. Assume that there exists x0 ¼ j2SB lj Aj a convex combination of the columns Aj ; j 2 S B which is compatible with Q ZB . As fAj ; j 2 S B g is minimal (Theorem 3.1), then lj – 0 for all j 2 S B . The columns xv  and x 0 are compatible, so that there exist kand P P P P P P k 0 such that j2SB v j Aj ¼ j2P B kj Aj and j2SB lj Aj ¼ j2P B k0 Aj . j Let j0 2 S B

X kj

Aj0 ¼

j2P B

¼

v j

In this section, we show that the objective function values of the reduced problems in IPS-A and IPS-N are equal at the first iteration. Let B be a feasible basis and RPB the reduced problem. Assume in this section that zSDB < 0 and let S B be the subset of the indices j 2 I B for which v j > 0. At the first iteration of IPS-N, we construct a new basis B0 . This basis is defined as B augmented by m  p* artificial variables combined with the p* variables xj ; j 2 P B . We add the incompatible variables xj ; j 2 S B to RP B0 . The new problem is denoted by RPDB . When we optimize this problem, we obtain a new objective function value zDRP . With the same reduced problem RPB and at the first iteration of IPS-A, we add the aggregate column P xv  ¼ j2SB v j Aj to the reduced problem. When we optimize this

j2S B nfj0 g

0

v j

Aj

0

X k0j X X lj Aj  Aj : l l j0 j2P  j0 j2S  nfj g B

B



0

From [7, Proposition 5.4], the columns fAj ; j 2 P B [ S B n fj0 gg lj

are linearly independent, so that k0j ¼ v 0 kj for j 2 P B and j0 P P lj   0 lj ¼ v  v j for j 2 S B . Since j2SB v j ¼ 1 and j2S B lj ¼ 1, we have j0

l = v*. A Theorem 5.2. zD RP ¼ zRP .

Proof. We prove that the optimal solution to RPDB is a feasible solution to RP AB and vice versa. Let xA be an optimal solution of RPAB and xD be an optimal solution to RPDB . We have xA ¼ ðxAC ; xx Þ and   xD ¼ ðxDC ; xv Þ with xAC and xDC in RjPB j ; xx in R and xv 2 RjSB j . The D constraint (2.5) of RP B is

X

xDj Aj þ

j2P B

X

xvj Aj ¼ b 

ð5:1Þ

j2S B

and the constraint (2.5) of RP AB is

X

xAj Aj þ xx xv  ¼ b;

ð5:2Þ

j2P B

P with xv  ¼ j2SB v j Aj . It is clear that ðxAC ; xx v  Þ is a feasible solution to (5.1), thus D zRP 6 zARP . We show now that there exists a positive k such that A ðxD C ; kÞ is a feasible solution of RP B . We have

X

xDj Aj þ

j2P B

where k ¼ 5. Relation between the iterations of IPS-A and IPS-N

X v j

X

Aj 

X

xvj Aj ¼ 

j2S B

P

X

xDj Aj þ kA0 ¼ b;

ð5:3Þ

j2P B

v  and A0 ¼ P j2S

j2S B xj

 xjv B

k

Aj . From (5.3) and since

Aj ; j 2 P B and b are compatible, then A0 is compatible. As   P xjv xjv  j2S B k ¼ 1 and from Lemma 5.1 we conclude that v j ¼ k , A0 ¼ xv  . Consequently ðxDC ; kÞ is a feasible solution of RPAB . Hence RP zRP h D P zA . RP In the previous theorem, we showed that zRP D ¼ zA at the first iteration. This equality is not preserved if we decide to add another convex combination x 2 XB , aggregated in RP AB and disaggregated RP in RP DB . We can only guarantee zRP D 6 zA in the second iteration.

Remark 5.3. Let x be an element of XB . In IPS-N, we construct a new basis B0 is Step 4 in order to add the columns forming x to RPD B . This means that we add some of the removed constraints to

42

A. Metrane et al. / European Journal of Operational Research 207 (2010) 37–44

Table 6.1 Instance characteristics. Instance

nb const

nb var

deg (%)

Instance

nb const

nb var

deg (%)

FA6 FA8 FA10 FA12 FA14 FA16 FA18 vcs1 vcs3 vcs5 vcs7 vcs9

5067 5159 5159 5159 5159 5182 5182 2084 2084 2084 2084 2084

17594 21437 24492 24645 22641 23990 24517 10343 6766 7837 8795 10150

68 65 66 66 71 64 65 44 45 48 47 50

FA7 FA9 FA11 FA13 FA15 FA17 FA19 vcs2 vcs4 vcs6 vcs8 vcs10

5159 5159 5159 5159 5182 5182 5182 2084 2084 2084 2084 2084

20434 23258 24812 25746 23650 24282 24875 6341 7337 8308 9241 6327

59 66 66 65 63 65 65 45 48 48 47 45

6. Computational experiments We have tested the IPS-A algorithm on degenerate linear programs and compared its performance to the IPS-N algorithm and the primal simplex of CPLEX solver, version 9.1.3, hereafter abbreviated as CPLEX. All tests were performed on a linux PC with a processor of 2.8 GHz. For the IPS-A and IPS-N algorithms, we used the same solver for the solution of the reduced problems RPB, and the dual simplex of CPLEX solver for the solution of the complementary problems SDB . For reasons of simplicity, the computation of the compatibility matrices was performed outside of CPLEX using the UMFACK library, version 5.0.3. 6.1. Instances The instances used in our computational experiments were derived from two types of problems, namely, the combined fleet assignment and aircraft routing (FA) problems and the simultaneous vehicle and crew scheduling (VCS) problem in public transit. These problems were selected because their linear relaxations are often highly degenerate when solved by column generation. The FA problem consists of assigning one aircraft type to each flight segment of a cyclic weekly schedule so as to maximize anticipated profits while building aircraft routes that respect maintenance requirements (each aircraft must be maintained at least every four days) and aircraft availability per type. The FA instances were generated from a real-world data set containing 2505 flight segments and four aircraft types. To do so, we used the column generation model of Barnhart et al. [1] which contains one variable for each valid sequence of flights (called strings) from the end of one maintenance check to the end of another one. In this model, there is one set partitioning constraint for each flight segment, one availability constraint for each aircraft type, and flow conservation constraints between the strings. Again, to obtain various FA instances, we applied a column generation algorithm to the available data set and collected different restricted master problems. These restricted master problems form our test set FA6 to FA19. Our instances have a higher proportion of degenerate variables than those of [6]. The VCS problem consists of simultaneously building driver workdays and bus schedules to cover a set of bus trip segments

at minimum cost. To generate VCS instances, we used the random generator of Haase et al. [9] and their column generation model which contains one set partitioning constraint for each trip segment and additional constraints to count the maximum number of buses used during a day. In this model, one variable is used to count the number of buses and there is one variable for each feasible driver workday generated by the column generator. To obtain different VCS instances, we ran a column generation algorithm on different random data sets and gathered various restricted master problems arising at different iterations of the column generation process. The number of constraints (nb const) and the number of variables (nb var) in each of these instances are given in Table 6.1. For each instance, we also specify the average percentage of degenerate basic variables (deg) in the reduced problem solutions computed throughout the IPS algorithm [6]. The higher this percentage is, the more seriously degeneracy impedes the solution process when the CPLEX solver is applied to solve an instance. 6.2. Mixed improved primal simplex It is well known (see Lasdon [10]) that in the column generation method the objective function decreases rapidly in the first iterations and slows downs as we are close to the optimal solution. This phenomenon is also observed when solving instances vcs4 (Fig. 1) and FA11 (Fig. 2) by IPS-A. To overcome this difficulty, we propose a hybrid algorithm referred to by IPS-M. It is a mixture of IPS-A and IPS-N. The key idea is to start with IPS-A and end up with IPS-N. Let  > 0 and let zRP be the objective value of the reduced problem at i iteration i. We keep generating the aggregate columns by solving RP SDB until jzRP i  ziþ1 j <  or the number of aggregate columns is lar-

8750 IPS-A IPS-N

8700

objective function value

RPD B to make these columns compatible with the new basis. In contrast to that, in IPS-A the aggregate column x is added to RP AB but without modifying the basis since x is compatible with it. Naturally, this implies that it takes longer to add a convex combination with negative reduced cost in IPS-N than in IPS-A, but the objective function decreases faster in IPS-N than in IPS-A if we add many elements of XB simultaneously to the reduced problem.

8650 8600 8550 8500 8450 8400 8350 10

15

20

25

30

35

40

45

50

time (s) Fig. 1. Comparison of the objective function decrease of IPS-A and IPS-N for vcs4.

43

A. Metrane et al. / European Journal of Operational Research 207 (2010) 37–44

-1.2945e+07 IPS-A IPS-N

objective function value

-1.295e+07 -1.2955e+07 -1.296e+07 -1.2965e+07 -1.297e+07 -1.2975e+07 -1.298e+07 -1.2985e+07 -1.299e+07 15

20

25

30

35

40

45

50

55

60

time (s) Fig. 2. Comparison of the objective function decrease of IPS-A and IPS-N for FA11.

ger than a certain threshold. This threshold is computed using the number of constraints of RPB. We then switch to IPS-N and continue on until the solution process is ended. A more sophisticated criterion may provide better results. 6.3. Results To compare the CPLEX, IPS-A, IPS-N and IPS-M algorithms, we conducted the following three-step experiment for each instance listed in Table 6.1. First, a feasible basis is obtained by executing

phase I of the primal simplex of CPLEX. Second, starting from this basis, phase II of the primal simplex of CPLEX is performed to reach optimality. Third, starting from the same feasible basis, the IPS-A, IPS-N and IPS-M algorithms are applied to compute an optimal solution. Because these algorithms must begin with a feasible basis, the time required to compute the feasible basis is not taken into account in any of the comparative results presented below including those concerning the primal simplex of CPLEX. Note that, for all tests, the default parameter setting of CPLEX was used when solving the different linear problems arising in these algorithms. The results of our experiments are reported in Table 6.2 for FA instances and Table 6.3 for VCS instances. For each instance, this table reports the total computational time (all times are in seconds) of CPLEX, IPS-N, IPS-A and IPS-M algorithms. The three last columns of this table provide the ratio of the CPLEX computational time over the IPS-N, IPS-A and IPS-N total computational time. After the results for each table, a row provides the averages of these statistics. Note that, for the IPS-N, IPS-A and IPS-M algorithms, the total computational time is the sum of the times devoted to solving the reduced and complementary problems and corresponds essentially to the time needed to compute the compatibility matrices (by UMFPACK) and set up the reduced and complementary problems. For the FA instances, we observe in Table 6.2 that IPS-N reduces the computational time by an average factor of 10.71 when compared to CPLEX. The IPS-A gained an average of 10.83 compared to CPLEX even if the number of iterations of IPS-N is smaller than for the IPS-A. This is not surprising because the time to add the incompatible columns that form the compatible convex combination with negative reduced cost x 2 XB in IPS-N is longer than

Table 6.2 Computational results for FA instances. inst

CPLEX

FA6 FA7 FA8 FA9 FA10 FA11 FA12 FA13 FA14 FA15 FA16 FA17 FA18 FA19 Avg

273.59 363.60 405.50 517.40 697.85 553.56 610.05 694.48 646.32 571.60 927.65 602.55 649.26 624.51

IPS-N

IPS-A

IPS-M

time

nb it

time

nb it

time

nb it

41.75 53.28 52.88 59.94 52.98 54.99 58.98 64.70 77.89 41.30 47.03 63.11 59.42 45.68

10 7 8 9 7 7 11 8 14 7 7 10 8 5 8.57

49.20 55.64 45.47 51 58.94 57.43 60.04 63.74 60.53 44.14 46 65.24 58.33 45.58

11 9 10 11 10 11 7 10 15 7 7 10 8 9 9.64

21.02 42.79 36.08 48.09 42.59 49.23 54.19 54.86 49.40 36.10 39.64 63.75 49.48 43.85

10 9 7 9 7 6 9 11 8 5 5 13 7 8 8.14

CPLEX IPSN

CPLEX IPSA

CPLEX IPSM

6.55 6.82 7.67 8.63 13.17 10.07 10.34 10.73 8.30 13.84 19.72 9.55 10.93 13.63 10.71

5.56 6.53 8.92 10.15 11.84 9.64 10.16 10.90 10.71 12.95 20.17 9.24 11.13 13.67 10.83

13.02 8.50 11.24 10.76 16.39 11.24 11.26 12.66 13.08 15.83 23.40 9.45 13.12 14.24 13.16

Table 6.3 Computational results for vcs instances. inst

CPLEX

vcs1 vcs2 vcs3 vcs4 vcs5 vcs6 vcs7 vcs8 vcs9 vcs10 Avg

235.47 97.71 109.08 134.25 168.86 164.57 189.88 206.60 247.97 102.26

IPS-N

IPS-A

IPS-M

time

nb it

time

nb it

time

nb it

63.11 36.73 36.85 40.17 44.27 45.74 49.13 49.45 60.94 31.12

10 8 8 10 11 9 11 10 10 7 9.4

73.22 38.83 41.99 49.29 55.79 57.44 59.07 62.17 76.16 38.64

13 11 11 13 15 14 15 14 16 11 13.4

66.34 34.80 40.55 43.23 49.62 52.57 52.22 50.21 66.75 33.21

11 9 9 12 11 11 11 10 11 8 10.3

CPLEX IPSN

CPLEX IPSA

CPLEX IPSM

3.73 2.66 2.96 3.34 3.81 3.60 3.86 4.18 4.07 3.29 3.55

3.22 2.52 2.60 2.72 3.03 2.87 3.21 3.32 3.26 2.65 2.94

3.55 2.81 2.69 3.11 3.40 3.13 3.64 4.11 3.71 3.08 3.32

44

A. Metrane et al. / European Journal of Operational Research 207 (2010) 37–44

adding x directly to IPS-A (see Remark 5.3). The IPS-M increases the average factor and decreases the number of iterations compared to IPS-N and IPS-A. The results show an impressive factor of 13.16 over CPLEX. This explains that the objective function value of IPS-A decreases rapidly compared to the objective function value of IPS-N for the first iterations but it decreases more slowly as soon as it approaches the optimal solution. For the VCS instances, we observe in Table 6.3 that IPS-N reduces the computational time by an average factor of 3.55 when compared to CPLEX. The IPS-A gained an average of 2.94 compared to CPLEX. This gain is less than in IPS-N because the number of iterations is much greater in IPS-A (13.4) than in IPS-N (9.4). The time per iteration of IPS-A is not reduced enough compared to IPS-N and IPS-N is faster than IPS-A for these problems. With IPS-N the number of iterations decreases compared to IPS-A but remains larger than with IPS-N. The average reduction factor for IPS-M is larger than in IPS-A but is smaller than in IPS-N. The less good performance of IPS-M for these problems can be explained by Figs. 1 and 2. The interval where IPS-A is in advance on IPS-N is smaller for VCS than for FA problems and the advance of IPS-A is also smaller. Furthermore, the objective function of IPS-N continues to decrease at a good rate after this interval. It is more difficult for VCS problems then for FA to switch from IPS-A to IPS-N at a good time and obtain improved results with the hybrid algorithm. 7. Conclusions In this paper, we present an interpretation of IPS algorithm recently introduced by Elhallaoui et al. [7] as a Dantzig–Wolfe decomposition. The master problem is the reduced problem of IPS containing mainly non degenerate constraints. The subproblem is the complementary problem containing degenerate constraints. It is a non aggregated decomposition in which the set of non zero variables of the optimal extreme point of the subproblem are added to the master problem at each iteration. This new interpretation of IPS as a Dantzig–Wolfe decomposition opens up the way to adapt to IPS the big amount of knowledge developed during the last 30 years to speed-up the Dantzig–Wolfe decomposition. This interpretation permits to introduce IPS-A, a variant of IPS corresponding to the classical Dantzig–Wolfe decomposition in which extreme points of the subproblem are added to the master problem. This new variant is faster than IPS on some problem and slower on others. This paper presents also an hybrid algorithm IPS-M combining the two variants of IPS-A and IPS-N. In average, this hybrid algorithm is up to 13.16 times faster than CPLEX on some types of problems. More sophisticated variants and strategies can be explored in the future to obtain a faster algorithm. Experiments show that IPS and its variants are efficient when the percentage of degeneracy is between 25% and 75%. In other

words, these algorithms are efficient when the reduced problem and the complementary problem are significantly smaller than the original problem. As a potential application, we would like to outline that in integer linear programming, often problems are modeled using many constraints in order to tighten the lower bound. Furthermore, cuts are usually added to the problem during the solving process for the same reason. Obviously, a large number of constraints slow down the solving process. Note that many of the constraints become degenerate. We think that using IPS or its variants will significantly reduce the solution time of integer linear problems. Finally, concepts developed in this paper have already been specialized for solving efficiently set partitioning problems. These problems are well known to be highly degenerate. This specialization allowed to reduce solution time by factors up to 45 for linear relaxation and up to 100 for integer solution (see [4–6]). Two possible research avenues may be the extension of these concepts to reduce degeneracy in the dual simplex and nonlinear programming algorithms.

References [1] C. Barnhart, N. Boland, L.W. Clarke, E.L. Johnson, G.L. Nemhauser, R.G. Shenoi, Flight string models for aircraft fleeting and routing, Transportation Science 32 (1998) 208–220. [2] G.B. Dantzig, Linear Programming and Extensions, Princeton University Press, Princeton, NJ, 1963. [3] G.B. Dantzig, P. Wolfe, Decomposition principle for linear programs, Operations Research 8 (1960) 101–111. [4] I. Elhallaoui, D. Villeneuve, F. Soumis, G. Desaulniers, Dynamic aggregation of set partitioning constraints in column generation, Operations Research 53 (2005) 632–645. [5] I. Elhallaoui, A. Metrane, F. Soumis, G. Desaulniers, Multi-phase dynamic constraint aggregation for set partitioning type problems, Mathematical Programming Series A (2008), doi:10.1007/s10107-008-0254-5. [6] I. Elhallaoui, G. Desaulniers, A. Metrane, F. Soumis, Bi-dynamic constraint aggregation and subproblem reduction, Computers & Operations Research 35 (2008) 1713–1724. [7] I. Elhallaoui, A. Metrane, G. Desaulniers, F. Soumis, An improved primal simplex algorithm for degenerate linear problems, Informs Journal on Computing, submitted for publication. [8] P.C. Gilmore, R.E. Gomory, A linear programming approach to the cutting stock problem–Part II, Operations Research 11 (1963) 863–888. [9] K. Haase, G. Desaulniers, J. Desrosiers, Simultaneous vehicle and crew scheduling in urban mass transit systems, Transportation Science 35 (2001) 286–303. [10] L.S. Lasdon, Optimization Theory for Large Systems, MacMillan, New York, 1970. [11] A. Löbel, Experiments with a Dantzig–Wolf decomposition for multiple-depot vehicle scheduling problems, ZIB, Preprint SC 97-16. [12] P.-Q. Pan, A basis deficiency-allowing variation of the simplex method for linear programming, Computers and Mathematics with Applications 36 (3) (1998) 33–53. [13] V. Raymond, F. Soumis, D. Orban, A new version for the IPS algorithm, Computers and Operations Research 37 (1) (2010) 91–98. [14] J.M. Valério de Carvalho, Exact solution of bin-packing problems using column generation and branch-and-bound, Annals of Operations Research 86 (1999) 629–659.