Systems Optimization: Column Generation

Column Generation Teo Chung-Piaw (NUS) 1 Lecture 1.1 Outline Slide 1 • Cutting Stock Problem • Classical Integer Programming Formulation • Set Covering Formulation • Column Generation Approach • Connection with Lagrangian Relaxation • Computational issues 2 Cutting Stock Problem 2.1 2.1.1 Introduction Example • A paper company has a supply of large rolls of paper, each of width W . Slide 2 Slide 3 • Customers demand ni rolls of width wi (i=1, . . . , m). (wi ≤ W ) Example: Quantity Ordered ni 97 610 395 211 Order Width (inches) wi 45 36 31 14 Slide 4 • The demand can be met by slicing a large roll in a certain way, called a pattern. • For example, a large roll of width 100 can be cut into 1 – 4 rolls each of width 25, or – 2 rolls each of width 35, with a waste of 30. 3 Solution Approach I 3.1 3.1.1 L. V. Kantorovich Formulation Slide 5 (1939 Russian, 1960 English) “Mathematical Methods of Planning and Organising Production” Management Science, 6, 366-422. • K: Set of available rolls. • y k : 1 if roll k is cut, 0 otherwise. • xki : number of times item i is cut on roll k. Objective: To minimize the number of rolls used to meet all the demand
min yk k∈K Slide 6 Constraints • Total number of times item i is cut is not less than the demand.
xki ≥ ni k∈K • The width of a roll is at most W
wi xki ≤ W y k i 2 min s.t. Slide 7  k∈K  k∈K  xki ≥ ni yk for a fixed item i, k i wi xi k ≤ W y for fixed roll k, xki ≥ 0, 0 ≤ y k ≤ 1 Integrality constraints on all variabes 3.1.2 Quality of Solution Slide 8 Scenario I: N1 • ni : uniform, between 1 and 100 (rand(100)+1); • wi : uniform, between 1 and 30 (rand(30)+1); • Width of Roll, W = 3000; Rolls 30 50 100 200 Items 60 100 200 400 constr 90 150 300 600 Note 1 variables 1830 5050 20100 80200 CPU (s) 2.8 14.33 179 3048 Code The OPL code for the problem: 3 Slide 9 Scenario II: Change width of the roll from 3000 to 150. Rolls Items constr variables CPU (s) 70 10 80 770 3.18 - never 140 20 160 2940 58.28 - never 210 30 240 6510 Out of memory The performance has deteriorated. Why? How good is the LP relaxation?  wi ni Observation: ZLP = iW . N2 This bound is trivial: The objective is to
min yk Slide 10 k∈K the optimal solution will satisfy: • Choose y k as small as possible. Therefore
wi xki = W y k i for all k Slide 11 • Choose xki as small as possible. Therefore
xki = ni k 4 for all i.
yk = k∈K =
 wi xk

wi i i = xki W W i k∈K k∈K
wi

wi ni k xi = W W i i k∈K Note 2 Proof We have the constraints
wi xki ≤ W y k . i  In the LP, since the objective is to minimize k y k , the optimal LP solution will be such that  k i wi xi = yk . W y k will be small if the xki values are small. At the same time, because
xki ≥ ni , k∈K to make xik small, at the optimal solution, we must have
xki = ni . k∈K So the objective function, for the optimal LP solution, reduces to

 wi xk i i yk = W k∈K k∈K

wi = xik W i k∈K
wi
= xik W i k∈K
wi ni = W i Slide 12 Another Example: W = 273 Quantity Ordered 233 310 122 157 120 Order Width (inches) 18 91 21 136 51 5 • LP: solved in 0.27s. Solution 228.7106. • IP: halted after 12 hours of computational time! 4 Approach II 4.1 4.1.1 Gilmore and Gomory Set Covering P. C. GILMORE AND R. E. GOMORY, A linear programming approach to the cutting-stock problem, Oper. Res., 8 (1961), pp. 849-859. xj = number of times pattern j is used aij = number of times item i is cut in pattern j For example, a large roll of width 100 can be cut into Slide 13 • 4 rolls each of width wi = 25 (pattern j, aij = 4) • 2 rolls each of width wk = 35 (pattern l, akl = 2) min s.t. Slide 14 n xj i = 1, . . . , m, j=1 aij xj ≥ ni , xj ∈ Z + , j = 1, . . . , n n j=1 Slide 15 Example: An instance: W = 100, m = 3. wi 25 35 45 1 4 0 0 2 2 1 0 pattern 3 4 2 1 0 2 1 0 5 0 1 1 6 0 0 2 ni 150 200 300 Slide 16 Another way to formulate the cutting stock problem: x1 4x1 + 0x1 + 0x1 + minimize x2 x3 x4 2x2 + 2x3 + 1x4 + 1x2 + 0x3 + 2x4 + 0x2 + 1x3 + 0x4 + 6 x5 0x5 + 1x5 + 1x5 + xj x6 0x6 ≥ 0x6 ≥ 2x6 ≥ j=1 RHS 150 200 300 xj = number of rolls to be cut using pattern j. 4.1.2 Computational Issues Slide 17 Linear relaxation: (LP) min s.t. n x n j=1 j j=1 aij xj ≥ ni , i = 1, . . . m, xj ≥ 0, j = 1, . . . , n 6 • LP solution provides a lower bound to IP. How to solve (LP)? Slide 18 • Feasible patterns: all nonnegative integer vectors (z1 , . . . , zm ) satisfying m
wi zi ≤ W i=1 • Not all feasible patterns are needed in the above formulation. Issue: In the constraint matrix A, the number of decision variables (i.e., number of feasible patterns) is large! 5 Cutting Stock Problem 5.1 5.1.1 Column Generation Algorithm 1. Start with a basic feasible solution B. Slide 19 For example, use the simple pattern to cut a roll into W/wi  rolls of width wi . (The basis matrix is a diagonal matrix.) m 2. For any pattern j, reduced cost is 1 − i=1 πi aij , where (π1 , . . . , πm ) (=cB B−1 ) is the simplex multipliers vector associated with the current basis. 3. Identify a pattern with negative reduced cost, or prove that none exists. Update basis and repeat. Slide 20 For instance, we may want to find the column with most negative reduced cost: m Z ∗ (π) = min 1 − i=1 πi xi  subject to i wi xi ≤ W, xi integral. • If Z ∗ (π) ≥ 0, all columns have negative reduced cost! • Otherwise, the solution gives rise to a column with negative reduced cost! 5.1.2 Identifying Columns m min subject to 1 − i=1 πi xi  i wi xi ≤ W, xi integral. is equivalent to solving  m Z  (π) =  max i=1 πi xi subject to i wi xi ≤ W , xi integral. 7 Slide 21 5.1.3 Knapsack problem Z � (π) = max subject to m πi x i  i=1 i Slide 22 wi xi ≤ W, xi integral. • The column generation method depends critically on how fast we can solve the knapsack problem. • How difficult is it to solve the knapssack problem? Computational Result on random instances using the MIP solver from CPLEX: n 1,000 10,000 100,000 Slide 23 CPU (s) 0.22 1.04 75.52 More specialized algorithm can be used to solve the Knapsack problem efficiently in practice. 5.1.4 How Good is the bound? Slide 24 Number of items = 5. Tested on several instances: Optimal 15 11 34 19 Col Gen (LP) 14.0533 10.4733 33.0989 18.2867 Round Up Conjecture: ZIP ≤ ZLP ? Slide 25 Unfortunately, this is not true: • W = 273 • w1 = 18 w2 = 91 w3 = 21 w4 = 136 w5 = 51 n1 n2 n3 n4 n5 = 233 = 310 = 122 = 157 = 120 8 ZLP (CG) = 228.9982. ZIP = 230. N3 Note 3 Round Up Conjecture In 1985 the theoretical result of Marcotte (The cutting stock problem and integer rounding, Math. Programming, 33 (1985), pp. 82-92) shed light on the relationship between solutions of the LP relaxation and the cutting stock problem itself. She proved that for some practical instances of the problem a so-called round-up property is valid. It means that to find the optimal value of the cutting stock problem, it is sufficient just to solve the LP relaxation and to round up the value of the objective function. Unfortunately, this conjecture is not true for all instances of the cutting stock problem. Fieldhouse (The duality gap in trim problems, SICUP-Bulletin No. 5, 1990) presents an example of the cutting stock problem with a gap of 1.0333. This gives rise to the “modified” round up conjecture. 5.1.5 Modified Round Up Conjecture ZIP ≤ ZLP  + 1? Slide 26 • This conjecture has not been answered. • Can you disprove it? 5.1.6 Getting Intergal Solution The solution obtained from solving the column generation problem may fractional. How to obtain integral solution? Slide 27 - round up the fractional solution. (e.g., change 18.3 to 19). - round down the fractional solution, and resolve the problem with smaller set of demand. - branch and bound to obtain the optimal integral solution 5.1.7 Rounding Up • Let xj be the (fractional) LP solution obtained from the column generation method. • Let xj = xj . xj integral.     • j ai,j xj ≥ ni implies j ai,j xj ≥ ni . So xj defined in this way is a feasible integral solution. • How good is this heuristic? 9 Slide 28 Slide 29 W = 273 w1 = 18 w2 = 91 w3 = 21 w4 = 136 w5 = 51 n1 n2 n3 n4 n5 = 233 = 310 = 122 = 157 = 120 • ZLP (CG) = 228.9982. • Round-Up produces a solution of 231 Slide 30 cut cut cut cut cut cut cut Round-Up 0 104 9 79 0 24 15 Fractional 0.0000 103.3333 8.2363 78.5000 0.0000 24.0000 14.9286 18 15 0 0 0 0 1 14 91 0 3 0 0 0 0 0 21 0 0 13 0 0 0 1 136 0 0 0 2 0 0 0 51 0 0 0 0 5 5 0 Disadvantages of the round-up heuristic? 6 6.1 6.1.1 Column Generation Dual Perspective Lagrangian Relaxation How would you use LR to solve the cutting stock problem? K k min k=1 y K s.t. xk ≥ ni ∀ i = 1, 2, . . . , m, mk=1k i k ∀ k = 1, . . . , K, i=1 xi wi ≤ W y Slide 31 y k ∈ {0, 1} ∀ k, xki ≥ 0, xki integral Which constraints would you relax? Suppose you relax the first class of constraints: L(u) = min s.t. K k   K k i=1 ui ni − k=1 xi m k=1 y + i=1 xki wi ≤ W y k ∀ k = 1, . . . , K, m y k ∈ {0, 1} ∀ k, xki ≥ 0, xki integral 10 Slide 32 where ui ≥ 0 for all i. Slide 33 L(u) = K
Lk (u) + u i ni i=1 k=1 where m
m y k − i=1 ui xik m k k s.t. i=1 xi wi ≤ W y k y ∈ {0, 1} Lk (u) = min xki ≥ 0, xki integral Lk (u) is the minimum of the two values: zero (when y k = 0), or 1-Z ∗ (when y k = 1) , where Z ∗ is obtained by solving the knapsack problem m Z ∗ = max ui xki m i=1k s.t. i=1 xi wi ≤ W Slide 34 xki ≥ 0, xki integral The subproblem in Lagragian Relaxation reduces again to a Knapsack problem! Slide 35 Proposition: N4 • maxu≥0 L(u) = ZLP (CG). • The optimalLanagrangian multiplers is the LP dual multiplers to the N constraints j=1 ai,j Kλj ≥ ni in the column formulation. • Lagrangian relaxation solves the dual of the column formulation! Note 4 Derivation of proposition L(u) = K
Lk (u) + m
ui ni i=1 k=1 The value Lk (u) does not depend on the index k. L(u) = KL(u) + m
ui ni i=1 where L(u) = min s.t. y− m m i=1 ui xi xi wi ≤ W y y ∈ {0, 1} i=1 xi ≥ 0, xi integral 11 • Let zj = (a1,j , . . . , am,j ), j = 1, . . . , N be the extreme points of the polytope m
xi wi ≤ W, xi ≥ 0, xi integral. i=1 • L(u) = min(0, 1−maxj=1,...,N m i=1 ai,j ui ) = minj=1,...,N min(0, 1− L(u) = KL(u) + m
m i=1 ui ni i=1 The Lagrangian Dual max L(u) = u≥0  m
max min ai,j ui ) K min(0, 1 − u≥0 j=1,...,N + m
 i=1 ui ni i=1 reduces to s.t. max y m y ≤ i=1 ui ni m m y ≤ K(1 − i=1 ui ai,j ) + i=1 ui ni for the zero extreme point for the jth extreme point ui ≥ 0 ∀ i = 1, . . . , m. Equivalently, max m y (λ0 ) y − i=1 ui ni ≤ 0 m m (λj ) y + K( i=1 ui ai,j ) − i=1 ui ni ≤ K ui ≥ 0 ∀ i = 1, . . . , m. λ0 , λj are the associated dual variables. The dual of this problem: N min j =1 Kλj N λ0 + j=1 λj = 1 N N −ni (λ0 + j=1 λj ) + j=1 ai,j Kλj ≥ 0. λj ≥ 0 ∀ j = 1, . . . , N. This is just min   j =1 Kλj N 12 ai,j ui ) λ0 + N N j=1 λj = 1 j=1 ai,j Kλj ≥ ni λj ≥ 0 ∀ j = 1, . . . , N. N For K large, the constraint λ0 + j=1 λj = 1 is redundant – as long as there N is a feasible solution λj with j =1 λj ≤ 1 Let xj = Kλj . We obtain the column formulation of the cutting stock problem ! 6.2 Lagrangian Relaxation 6.2.1 Comparison Slide 36 • The bounds obtained by both methods are identical. • Which method is better? Column Generation Primal, dual optimal solution (Primal) Bounds monotone Dual solution zig-zag LP solver needed 6.3 Langrangean Relaxation Dual but not primal solution (Dual) Bounds zig-zag Dual solution suitably selected Easy to implement Speeding Up 6.3.1 Column Selection Slide 37 • Prevent generation of redundant columns.   T Instead of minj cj − π Nj , solve min j:π T Nj >0 cj , π T Nj   cj − π T N j min . j 1 T Nj or Slide 38     c −π T N • minj cj − π T Nj versus minj j 1T Nj j cut cut cut cut cut cut cut Round-Up 0 104 9 79 0 24 15 Fractional 0.0000 103.3333 8.2363 78.5000 0.0000 24.0000 14.9286 18 15 0 0 0 0 1 14 13 91 0 3 0 0 0 0 0 21 0 0 13 0 0 0 1 136 0 0 0 2 0 0 0 51 0 0 0 0 5 5 0 • Maintain a column pool. Check for columns with negative reduced cost in the column pool before solving the pricing subproblem. Replenish the column pool everytime you solve a pricing subproblem. 6.3.2 Dual Selection Slide 39 Restrict domain of dual prices Use properties of optimal dual prices to restrict the domain. In the cutting stock problem, suppose the orders are ranked such that w1 < w2 < . . . < wm , then it is easy to see that the dual prices satisfy: π1 ≤ π2 ≤ . . . ≤ πm . Translating into the primal, this is equivalent to adding the following columns with zero cost to the primal problem:       1 0 0  −1   1         ···   0   −1   ···         ··· ,  0 , ... ···         ···   ···   0         ···   ···   1  0 0 −1 7 Column Generation 7.1 7.1.1 Applications Examples Other application of column generation: Slide 40 Vehicle Routing with time window or other types of constraint: • A set of m customers. • Each customer must be served within certain time window. • Find a set of routes to serve all the customers, so that each customer will be visited by a vehicle within the stipulated time window. Here each column may represent a feasible trip. One likely objective is to find minimal number of vehicles to cover all the customers. 14 Slide 41 Column generation phase now reduces to the following: Given a profit πi for each demand point, find a route that satisfies: Slide 42 • feasibility constraints (meet the time window constraints); • total profits accrued by serving the demand points on the route is maximum - a variant of the TSP problem. N5 Note 5 Papers • J. Desrosiers, Y. Dumas,F. Soumis & M. Solomon. Time Constrained Routing and Scheduling, Handbooks in OR & MS, 8 (1995) • G. Desaulniers et al. A Unified Framework for Deterministic Vehicle Routing and Crew Scheduling Problems T. Crainic & G. Laporte (eds) Fleet Management & Logistics (1998). 8 Conclusions • Column Generation has been successfully used to solve many large scale integer programming problem arising in the industry. • Able to handle large scale model that standard commercial MIP solver cannot handle. • Ability to solve the pricing subproblem efficiently is key to the approach • Connection between Column generation and Lagragian Relaxation • Non-linearities occuring in practical problems taken care of in the subproblem (next lecture) 15 Slide 43