Description: Shows how to implement column generation in JuMP for solving the cutting stock problem.

Author: Chiwei Yan

License: ## Cutting Stock Problem Using Julia/JuMP¶

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.

### Mathematical Formulation¶

Sets

• $I$ = set of order widths
• $J$ = set of patterns

Parameters

• $a_{ij}$ = number of rolls of width $i$ cut in pattern $j$
• $b_i$ = demand for order width $i$

Decision Variables

• $x_j$ = number of rolls cut using pattern $j$

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}

### Let's see how to formulate this problem in JuMP¶

In :
#import necessary packages and define model
using JuMP
using MathOptInterface
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

Out:
\begin{alignat*}{1}\min\quad & 0\\ \text{Subject to} \quad\end{alignat*}

We are now going to initialize a "restricted master problem" with only two variables, corresponding two cutting patterns:

• width (14,31,36,45), quantity (1,1,0,1), denoted as $x_1$
• width (14,31,36,45), quantity (0,0,2,0), denoted as $x_2$

[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}
In :
#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

Out:
\begin{alignat*}{1}\min\quad & x_{1} + x_{2}\\ \text{Subject to} \quad & x_{1} \geq 211\\ & x_{1} \geq 395\\ & 2 x_{2} \geq 610\\ & x_{1} \geq 97\\ & x_{i} \geq 0 \quad\forall i \in \{1,2\}\\ \end{alignat*}
In :
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.

In :
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".

In :
#model before adding new column
master

Out:
\begin{alignat*}{1}\min\quad & x_{1} + x_{2}\\ \text{Subject to} \quad & x_{1} \geq 211\\ & x_{1} \geq 395\\ & 2 x_{2} \geq 610\\ & x_{1} \geq 97\\ & x_{i} \geq 0 \quad\forall i \in \{1,2\}\\ \end{alignat*}

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:

• What's the coefficient for this new variable in the objective function?
• Which constraint does this new variable appear? With what coefficient?
In :
#column-wise adding new variable z
@variable(master, z>=0, objective=1, inconstraints=myCons, coefficients=pattern)

#look at the master problem again
master

Out:
\begin{alignat*}{1}\min\quad & x_{1} + x_{2} + z\\ \text{Subject to} \quad & x_{1} \geq 211\\ & x_{1} + 2.9999999999999996 z \geq 395\\ & 2 x_{2} \geq 610\\ & x_{1} \geq 97\\ & x_{i} \geq 0 \quad\forall i \in \{1,2\}\\ & z \geq 0\\ \end{alignat*}
In :
#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

### Now it's time to put all pieces together¶

In :
#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]

@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?

### How to obtain INTEGER solution?¶

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.

#### Method 1: Rounding¶

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:

In :
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?

#### Method 2: Branch-and-Bound on Root Node¶

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.

In :
#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?