**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**

- $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}$$In [1]:

```
#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
```

Out[1]:

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?

In [2]:

```
#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[2]:

In [3]:

```
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
```

**[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 [4]:

```
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))
```

**[Result Analysis]**

The reduced cost of this variable is $(1-3)=-2<0$. Add this new variable to the ** "restricted master problem"**.

In [5]:

```
#model before adding new column
master
```

Out[5]:

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 [6]:

```
#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
```

Out[6]:

In [7]:

```
#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
```

**[Result Analysis]**

We see that after adding a new variable, the objective value is reduced to 577.3

In [8]:

```
#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
```

[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:

In [9]:

```
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)
```

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.

In [10]:

```
#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)
```

We save one paper roll by using method 2

[Question]:Is method 2 always able to produce feasible integer solution?