Description: Shows how to implement column generation in JuMP for solving the cutting stock problem.
Author: Chiwei Yan
License:
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
In this notebook, we deploy column generation to solve the cutting stock problem using Julia/JuMP. The origin of the cutting stock problem is in the paper industry. Given paper rolls of fixed width and a set of orders for rolls of smaller widths, the objective of the cutting stock problem is to determine how to cut the rolls into smaller widths to fulfill the orders in such a way as to minimize the amount of scrap. The cutting stock problem example we use in this notebook is from Linear Programming by Vasek Chvatal, 1983.
Suppose that rolls are produced in a uniform width of 100 inches and that orders can be placed for rolls of widths 14 inches, 31 inches, 36 inches, and 45 inches. The company has received the following orders:
\begin{array}{|c|c|} \hline \textrm{Order Width} & \textrm{Quantity} \\\hline 14 & 211 \\\hline 31 & 395 \\\hline 36 & 610 \\\hline 45 & 97 \\\hline \end{array}A single 100 inch roll can be cut into one or more of the order widths. For example, one roll could be cut into two rolls of 45 inches with a 10 inch roll of scrap.
Or a roll could be cut into a roll of 45 inches, a roll of 31 inches, and a roll of 14 inches with no scrap. Each such possible combination is called a pattern. For this example, there are 37 different patterns. Determining how many of each pattern to cut to satisfy the customer orders while minimizing the scrap is too difficult to do by hand. Instead the problem can be formulated as an optimization problem, specifically an integer linear program.
Sets
Parameters
Decision Variables
The objective of the cutting stock problem is to minimize the number of rolls cut subject to cutting enough rolls to satisfy the customer orders. Using the notation above, the problem can be formulated as follows:
$$\begin{align} \nonumber \min\qquad &\sum_{j\in J}x_j \\ \nonumber \textrm{s.t.}\qquad &\sum_{j\in J}a_{ij}x_j\ge b_i,~\forall i\in I \\ \nonumber &x_j\in\textrm{integer},~\forall j\in J \end{align}$$#import necessary packages and define model
# Load JuMP
using JuMP
using MathOptInterface
# Load solver package
using MathOptInterfaceXpress
# shortcuts
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
master = Model(optimizer = XpressOptimizer())
#If Gurobi is installed (requires license), you may uncomment the code below to switch solvers
#using Gurobi
#master = Model(solver=GurobiSolver(Method=0)) # Switch LP algorithm to Primal Simplex, in order to enjoy warm start
We are now going to initialize a *"restricted master problem"* with only two variables, corresponding two cutting patterns:
[Recall 1] what's the meaning of each variable? Number of paper rolls cut using this pattern.
[Recall 2] How should the formulation of the restricted master problem look like?
$$ \begin{align} \nonumber\min\qquad\qquad\quad &x_1+x_2 \\ s.t.\qquad\left( \begin{array}{c} 1 \\ 1 \\ 0 \\ 1 \end{array} \right)&x_1+\left( \begin{array}{c} 0 \\ 0 \\ 2 \\ 0 \end{array} \right)x_2 \ge \left( \begin{array}{c} 211 \\ 395 \\ 610 \\ 97 \end{array} \right) \\ &x_1,x_2\ge0 \end{align}$$#define initial variables
@variable(master, x[1:2] >= 0)
#width
w=[14 31 36 45]
#constraint coefficient for initial variables
A=[1 0; 1 0; 0 2; 1 0]
b=[211; 395; 610; 97]
#define constraint references
myCons = Any[]
#define constraints
for i=1:4
myCon = @constraint(master, dot(x, vec(A[i,:]))>=b[i])
push!(myCons, myCon)
end
#define objective
@objective(master, Min, sum(x))
master
JuMP.optimize(master)
status = JuMP.terminationstatus(master)
JuMP.resultvalue.(x)
#get the optimal solution
println("\nOptimal Solution is:\n")
println("width: ", w)
epsilon=1e-6
for i=1:size(A,2)
if JuMP.resultvalue(x[i])>epsilon
println("Cutting Pattern: ", A[:,i], ", Number of Paper Rolls Cut Using this Pattern: ", JuMP.resultvalue(x[i]))
end
end
Optimal Solution is: width: [14 31 36 45] Cutting Pattern: [1,1,0,1], Number of Paper Rolls Cut Using this Pattern: 395.0 Cutting Pattern: [0,0,2,0], Number of Paper Rolls Cut Using this Pattern: 305.0
[Result Analysis]
The minimal number of paper rolls is 700.
Clearly this is not the best we can do, because we are not considering all possible feasible patterns.
Let's now generate some new patterns based on the value of reduced costs. Denote $r=(r_1,r_2,r_3,r_4)$ as the optimal dual price of constraints 1, 2, 3, 4. The reduced cost of a potential variable $x_k$, with cutting pattern $A_k$ can be calculated as $$rc(x_k)=1-A_k^Tr$$
We want to add a potential variable $x_k$ such that $rc(x_k)<0$, this can be done by solving the following sub-problem:
$$\begin{align} z^*=\max\qquad &r_1a_{k,1}+r_2a_{k,2}+r_3a_{k,3}+r_4a_{k,4} \\ s.t.\qquad &14a_{k,1}+31a_{k,2}+36a_{k,3}+45a_{k,4}\le 100 \\ &a_{k,1},a_{k,2},a_{k,3},a_{k,4}\ge0,~\textrm{and are integers} \end{align}$$If $z^*>1$, then $x_k$ with cutting pattern $(a_{k,1},a_{k,2},a_{k,3},a_{k,4})$ should be added to the master problem. And resolve the master problem.
r=[getDual(myCons)[1:4]]
sub = Model(optimizer = XpressOptimizer())
#width
w=[14,31,36,45]
#define cutting pattern variables
@variable(sub, a[1:4]>=0, Int)
#define feasible cutting constraint
@constraint(sub, dot(w,a)<=100)
#define objective
@objective(sub, Max, dot(r,a))
sub
JuMP.optimize(master)
status = JuMP.terminationstatus(master)
#print new cutting pattern
pattern=[getValue(a)[1:4]]
println("width: ", w)
println("\nNew Cutting Pattern: ", int(pattern))
width: [14,31,36,45] New Cutting Pattern: [0,3,0,0]
[Result Analysis]
The reduced cost of this variable is $(1-3)=-2<0$. Add this new variable to the *"restricted master problem"*.
#model before adding new column
master
JuMP supports column-wise modeling in defining variables. Think about it, when we add a new variable to the existing model, we need to know:
#column-wise adding new variable z
# TODO column addtion
@variable(master, z>=0, objective=1, inconstraints=myCons, coefficients=pattern)
#look at the master problem again
master
#solve the master problem again
JuMP.optimize(master)
status = JuMP.terminationstatus(master)
#get the optimal solution
println("\nOptimal Solution is:\n")
println("width: ", w)
for i=1:length(x)
if JuMP.resultvalue(x[i])>epsilon
println("Cutting Pattern: ", A[:,i], ", Number of Paper Rolls Cut Using this Pattern: ", JuMP.resultvalue(x[i]))
end
end
if JuMP.resultvalue(z)>epsilon
println("Cutting Pattern: ", int(pattern), ", Number of Paper Rolls Cut Using this Pattern: ", JuMP.resultvalue(z))
end
Optimal Solution is: width: [14,31,36,45] Cutting Pattern: [1,1,0,1], Number of Paper Rolls Cut Using this Pattern: 211.0 Cutting Pattern: [0,0,2,0], Number of Paper Rolls Cut Using this Pattern: 305.0 Cutting Pattern: [0,3,0,0], Number of Paper Rolls Cut Using this Pattern: 61.33333333333334
[Result Analysis]
We see that after adding a new variable, the objective value is reduced to 577.3
#import necessary packages and define master problem
using JuMP, Cbc
master = Model()
#If Gurobi is installed (requires license), you may uncomment the code below to switch solvers
#using Gurobi
#master = Model(solver=GurobiSolver(Method=0)) # Switch LP algorithm to Primal Simplex, in order to enjoy warm start
#define initial variables
@variable(master, x[1:2] >= 0)
#constraint coefficient for initial variables
A=[1 0; 1 0; 0 2; 1 0]
b=[211; 395; 610; 97]
#define constraint references (why?)
@constraint myCons[1:4]
#define constraints
for i=1:4
myCons[i] = @constraint(master, dot(x, vec(A[i,:]))>=b[i])
end
#define objective
@objective(master, Min, sum(x))
#solve master problem
JuMP.optimize(master)
println("Iteration 1, Master Problem Objective Value:", getObjectiveValue(master))
#subproblem to iteratively generate new columns
#get optimal dual prices from the master problem
r=[JuMP.resultdual(myCons)[1:4]]
sub=Model(optimizer = XpressOptimizer())
#width
w=[14,31,36,45]
#define cutting pattern variables
@variable(sub, a[1:4]>=0, Int)
#define feasible cutting constraint
@constraint(sub, dot(w,a)<=100)
#define objective
@objective(sub, Max, dot(r,a))
#solve the subproblem
JuMP.optimize(sub)
sub_obj=JuMP.objectivevalue(sub);
epsilon=1e-6;
#list of new variables
newColumns=Variable[]
#pattern list
A_new=Float64[];
iter=2
while sub_obj>1+epsilon #why?
#cutting pattern (constraint coefficients) for the new variable
pattern=JuMP.resultvalue(a)[1:4]
#column-wise adding new variable z
@variable(master, z>=0, objective=1, inconstraints=myCons, coefficients=pattern)
println("\tAdd a new variable with cutting pattern: ", pattern, ", reduced cost: ", (1-sub_obj))
#add new variable to the new variable list
push!(newColumns, z)
#add new cutting pattern to pattern list
append!(A_new, pattern)
JuMP.optimize(master)
println("\nIteration ",iter, ", Master Problem Objective Value:", JuMP.objectivevalue(master))
#get new optimal dual prices
r=[JuMP.resultdual(myCons)[1:4]]
#modify the objective of the subproblem based on new dual prices
@objective(sub, Max, dot(r,a))
JuMP.optimize(sub)
sub_obj=JuMP.objectivevalue(sub)
iter=iter+1
end
#print optimal solution
A_new=reshape(A_new,4, convert(Int64,length(A_new)/4))
println("\nOptimal Solution is:\n")
println("width: ", w)
for i=1:length(x)
if JuMP.resultvalue(x[i])>epsilon
println("Cutting Pattern: ", A[:,i], ", Number of Paper Rolls Cut Using this Pattern: ", JuMP.resultvalue(x[i]))
end
end
for i=1:length(newColumns)
if getValue(newColumns[i])>epsilon
println("Cutting Pattern: ", int(A_new[:,i]), ", Number of Paper Rolls Cut Using this Pattern: ", JuMP.resultvalue(newColumns[i]))
end
end
Iteration 1, Master Problem Objective Value:700.0 Add a new variable with cutting pattern: [0.0,2.9999999999999996,0.0,0.0], reduced cost: -2.0 Iteration 2, Master Problem Objective Value:577.3333333333334 Add a new variable with cutting pattern: [7.000000000000001,0.0,0.0,0.0], reduced cost: -3.666666666666666 Iteration 3, Master Problem Objective Value:517.6190476190477 Presolve 0 (-1) rows, 0 (-4) columns and 0 (-4) elements Optimal - objective value 1 After Postsolve, objective 1, infeasibilities - dual 0 (0), primal 0 (0) Optimal objective 1 - 0 iterations time 0.002, Presolve 0.00 Cbc0045I Solution with objective value -1 saved Add a new variable with cutting pattern: [2.0,0.0,2.0,0.0], reduced cost: -0.2857142857142856 Iteration 4, Master Problem Objective Value:501.33333333333337 Presolve 0 (-1) rows, 0 (-4) columns and 0 (-4) elements Optimal - objective value 1 After Postsolve, objective 1, infeasibilities - dual 0 (0), primal 0 (0) Optimal objective 1 - 0 iterations time 0.002, Presolve 0.00 Cbc0045I Solution with objective value -1 saved Add a new variable with cutting pattern: [0.0,0.0,0.0,2.0], reduced cost: -0.33333333333333326 Iteration 5, Master Problem Objective Value:485.1666666666667 Presolve 0 (-1) rows, 0 (-4) columns and 0 (-4) elements Optimal - objective value 1 After Postsolve, objective 1, infeasibilities - dual 0 (0), primal 0 (0) Optimal objective 1 - 0 iterations time 0.002, Presolve 0.00 Cbc0045I Solution with objective value -1 saved Presolve 0 (-1) rows, 0 (-4) columns and 0 (-4) elements Optimal - objective value 1 After Postsolve, objective 1, infeasibilities - dual 0 (0), primal 0 (0) Optimal objective 1 - 0 iterations time 0.002, Presolve 0.00 Cbc0045I Solution with objective value -1 saved Add a new variable with cutting pattern: [0.0,2.0,1.0,0.0], reduced cost: -0.16666666666666674 Iteration 6, Master Problem Objective Value:452.25 Optimal Solution is: width: [14,31,36,45] Cutting Pattern: [2,0,2,0], Number of Paper Rolls Cut Using this Pattern: 206.25 Cutting Pattern: [0,0,0,2], Number of Paper Rolls Cut Using this Pattern: 48.5 Cutting Pattern: [0,2,1,0], Number of Paper Rolls Cut Using this Pattern: 197.5 Presolve 0 (-1) rows, 0 (-4) columns and 0 (-4) elements Optimal - objective value 1 After Postsolve, objective 1, infeasibilities - dual 0 (0), primal 0 (0) Optimal objective 1 - 0 iterations time 0.002, Presolve 0.00 Cbc0045I Solution with objective value -1 saved Coin0505I Presolved problem not optimal, resolve after postsolve Coin0505I Presolved problem not optimal, resolve after postsolve
[Exercise and Discussion]:
- Change the initial variables we use to construct the first restricted master problem (but still maintain the starting restricted master problem feasible). How does it effect the convergence of the algorithm? (number of total columns generated?)
- Could you find a way to generate multiple columns whose reduced cost are less than 0 at one iteration?
We solve the LP relaxation successfully using column generation, however the original cutting stock problem is an integer program. Can we apply column generation to obtain optimal integer solution?
The answer is Yes. However, it involves an advanced solution methodology called branch-and-price where column generation is applied on each node of the branch-and-bound tree. Unfortunately, commercial solvers (Gurobi, CPLEX) don't support this feature. Till now, the only academic solver supports branch-and-price is SCIP.
Instead of solving the integer program to optimality, we here introduce two approximation methods that are widely used in solving real-world problems.
Rounding a fractional solution to its nearest and feasible is a common heuristic for solving integer program. It's pretty problem specific. In cutting stock problem, we observe that if we round up all the fractional solutions, feasibility will maintain. Thus we get our first integer solution:
println("\nInteger Solution Based on Rounding is:\n")
println("width: ", w)
summation=0
for i=1:length(x)
if JuMP.resultvalue(x[i])>epsilon
println("Cutting Pattern: ", A[:,i], ", Number of Paper Rolls Cut Using this Pattern: ", ceil(JuMP.resultvalue(x[i])))
summation=summation+ceil(JuMP.resultvalue(x[i]))
end
end
for i=1:length(newColumns)
if getValue(newColumns[i])>epsilon
println("Cutting Pattern: ", int(A_new[:,i]), ", Number of Paper Rolls Cut Using this Pattern: ", ceil(JuMP.resultvalue(newColumns[i])))
summation=summation+ceil(JuMP.resultvalue(newColumns[i]))
end
end
println("Total Number of Paper Rolls Used: ", summation)
Integer Solution Based on Rounding is: width: [14,31,36,45] Cutting Pattern: [2,0,2,0], Number of Paper Rolls Cut Using this Pattern: 207.0 Cutting Pattern: [0,0,0,2], Number of Paper Rolls Cut Using this Pattern: 49.0 Cutting Pattern: [0,2,1,0], Number of Paper Rolls Cut Using this Pattern: 198.0 Total Number of Paper Rolls Used: 454.0
We now have an integer solution using 454.0 paper rolls in total. Can we do better?
It is troublesome to implement column generation on every node of the branch and bound tree. A common industry / research practice is to directly branch-and-bound the model only with columns generated from solving the LP relaxation. This is a heuristic because optimal set of cutting patterns for the IP might not be the same as the LP relaxation, i.e. we might lose some "good columns" to reach optimal integer solution. The upside is, it is very easy to implement with commercial solvers.
#change the solver from Clp to Cbc in order to get support for integer variables
#if you use Gurobi as your solver choice, you don't need to switch solver.
setSolver(master, CbcSolver())
#change variable type from continuous to integer
for i=1:length(x)
setCategory(x[i], :Int)
end
for i=1:length(newColumns)
setCategory(newColumns[i],:Int)
end
JuMP.optimize(master)
#print optimal solution
println("\nInteger Solution Based on Branch-and-Bound is:\n")
println("width: ", w)
summation=0
for i=1:length(x)
if JuMP.resultvalue(x[i])>epsilon
println("Cutting Pattern: ", A[:,i], ", Number of Paper Rolls Cut Using this Pattern: ", JuMP.resultvalue(x[i]))
summation=summation+JuMP.resultvalue(x[i])
end
end
for i=1:length(newColumns)
if JuMP.resultvalue(newColumns[i])>epsilon
println("Cutting Pattern: ", int(A_new[:,i]), ", Number of Paper Rolls Cut Using this Pattern: ", JuMP.resultvalue(newColumns[i]))
summation=summation+JuMP.resultvalue(newColumns[i])
end
end
println("Total Number of Paper Rolls Used: ", summation)
Integer Solution Based on Branch-and-Bound is: width: [14,31,36,45] Cutting Pattern: [0,0,2,0], Number of Paper Rolls Cut Using this Pattern: 100.99999999999999 Cutting Pattern: [2,0,2,0], Number of Paper Rolls Cut Using this Pattern: 106.0 Cutting Pattern: [0,0,0,2], Number of Paper Rolls Cut Using this Pattern: 49.0 Cutting Pattern: [0,2,1,0], Number of Paper Rolls Cut Using this Pattern: 198.0 Total Number of Paper Rolls Used: 454.0 Presolve 0 (-4) rows, 0 (-7) columns and 0 (-11) elements Optimal - objective value 453 After Postsolve, objective 453, infeasibilities - dual 0 (0), primal 0 (0) Optimal objective 453 - 0 iterations time 0.002, Presolve 0.00 Cbc0045I Warning 3 integer variables were more than 1.0e-4 away from integer Cbc0045I Given objective value 452.25, computed 453 Cbc0045I Solution with objective value 453 saved
We save one paper roll by using method 2
[Question]:
Is method 2 always able to produce feasible integer solution?