## Cutting Stock Problem: A tutorial with example to understand Column generation

Column generation is an efficient technique to solve large-sized mathematical programs. It is also called delayed-column generation due to its nature to solve the problem. It is firstly studied and applied for the Cutting stock Problems, recognized by a Soviet mathematician and economist Leonid Kantorovich in 1939.

### Cutting Stock problem example

A Cutting stock problem is the problem of cutting off a minimum number of large-sized items to fulfill the customer requirements of various small-sized items of different quantities. For example, in the steel manufacturing industry, we have a stock of large-sized rods, each with a size of feet. The customer requires three different types of rods, 1) rods, each with the size of 4 feet, 2) rods of size feet, and 3) rods, each with the size of feet. Let us call them demand order types ,

The cutting stock problems are interesting – for the customer’s requirement with different items, problems become difficult. Suppose a customer asks for

cuts of; waste = cuts of; waste = cuts of; waste = cuts of and cut of; waste = cut of and cut of; waste =

And many such combinations of possibilities to meet the total demand.

## Idea -1: A dummy heuristic for Cutting Stock Problem

To solve this problem, the first (dummy) idea that comes to mind is to pick the large-sized rod and start cutting it until we produce the desired demand orders. The demand requirement is fulfilled one by one, i.e., first, we meet the

Let us do a simple python implementation for this –

```
large_size_rod_used = 0
demand_size = [4, 6, 7]
demand_count = [80, 50, 100]
large_len = 15
effective_len = large_len
waste = 0
for i in range(3):
while (True):
if (demand_count[i] > 0) :
if (demand_size[i] <= effective_len):
effective_len = effective_len - demand_size[i]
demand_count[i] = demand_count[i] - 1
else:
large_size_rod_used = large_size_rod_used + 1
waste = waste + effective_len
effective_len = large_len
continue
else:
break
print("large_size_rod_used ", large_size_rod_used)
print("total waste ", waste)
```

Let us run this:

We see that, to meet the demand, the total number of large-sized rods used is = 101, and waste accumulated in feet = 202 feet.

## Idea -II: : A dummy heuristic for Cutting Stock Problem

We can improve it by not just fulfilling the demand order by order. Instead, we can do the following: Start with the order type of minimum length so that the leftover part from the large rod can fulfill the other item type, which leads to minimizing waste.

Its python implementation:

```
import numpy as np
large_size_rod_used = 0
demand_size = [4, 6, 7]
number_of_demand_type = 3
demand_count = [80, 50, 100]
large_len = 15
effective_len = large_len
waste = 0
def indices_of_sorted_list(s):
return sorted(range(len(s)), key=lambda k: s[k])
# np_demand_size = np.array(demand_size)
for i in indices_of_sorted_list(demand_size):
while (True):
if (demand_count[i] > 0) :
if (demand_size[i] <= effective_len):
effective_len = effective_len - demand_size[i]
demand_count[i] = demand_count[i] - 1
while(True):
temp_remaining = [(effective_len - demand_size[j],j) for j in range(number_of_demand_type)
if effective_len >= demand_size[j] and demand_count[j]>0]
if temp_remaining:
j = np.argmin(np.array([x[0] for x in temp_remaining]))
j = [x[1] for x in temp_remaining][j]
effective_len = effective_len - demand_size[j]
demand_count[j] = demand_count[j] - 1
else:
break
large_size_rod_used = large_size_rod_used + 1
waste = waste + effective_len
effective_len = large_len
continue
else:
break
print("large_size_rod_used ", large_size_rod_used)
print("total waste ", waste)
```

Let us run this:

To meet the demand, the total number of large-sized rods used is = 95, and waste accumulated in feet = 105 feet. May be we are lucky and obtain the optimal solution for this particular problem. If we change the demand orders, this idea might not work. We can come up with such demand picture.

Since both methods do not use any technique that minimizes the number of large-sized rods, there is a high likelihood we will land up to the suboptimal solution.

## A basic mathematical modeling approach: a foundation of Column Generation for Cuttings Stock Problem

If we define

Similarly, suppose

Now, the constraints that can be obtained from the customer demand picture are

are binary, are positive integers.

The objective here is to minimize the sum of all the

Note that the estimated value of

If we consider

```
import cplex
import sys
# Input data
DEBUG = True
class DATA:
n= 100 # number of large-sized rod
m = 3 # number of required rod types
small_rod_demands = [80,50,100]
small_rod_size = [4,6,7]
large_rod_size = [15]
class CPLEX:
def __init__(self):
self.model = cplex.Cplex()
def solve(self):
self.model.solve()
def read_problem(self, problem):
return cplex.Cplex(problem)
def write(self, file_name):
self.model.write(file_name)
def get_objective_value(self):
return self.model.solution.get_objective_value()
def define_variable(self, varname_list, type_of_var, objective_coef_vector):
n_var = len(varname_list)
if type_of_var == "binary":
self.model.variables.add(obj=objective_coef_vector , names =varname_list,
types = [self.model.variables.type.binary] * n_var )
elif type_of_var == "continuous":
self.model.variables.add(obj= objective_coef_vector, names =varname_list,
types = [self.model.variables.type.continuous] * n_var)
elif type_of_var == "integer":
self.model.variables.add(obj= objective_coef_vector, names =varname_list,
types = [self.model.variables.type.integer] * n_var)
return None
def define_constraint(self, the_const_names, the_rows, the_senses, my_rhs ):
self.model.linear_constraints.add(lin_expr=the_rows, senses= the_senses, rhs = my_rhs, names = the_const_names)
return None
class cutting_stock:
def __init__(self, n, m, small_rod_demands,small_rod_size,large_rod_size):
self.n = n
self.m = m
self.small_rod_demands = small_rod_demands
self.small_rod_size = small_rod_size
self.large_rod_size = large_rod_size
self.bin_var_names = ["y_"+str(i) for i in range(n)] # binary variables
self.int_var_name_types = [["x_"+str(i)+"_"+str(j) for j in range(n)] for i in range(m)]
self.const_names_types = ["t_"+str(j) for j in range(m)]
self.const_sense_types = ["G" for j in range(m)]
self.const_rhs_types = small_rod_demands
self.const_name_pattern = ["p_"+str(i) for i in range(n)]
self.const_sense_pattern = ["L" for i in range(n)]
self.const_rhs_pattern = [0 for i in range(n)]
self.cplexSolver = CPLEX()
def create_mathematical_model(self):
# 1. binary variable
objective_coef_vector = [1 for i in range(self.n)]
self.cplexSolver.define_variable(self.bin_var_names, "binary", objective_coef_vector)
# 2. integer variable
for i in range(self.m):
objective_coef_vector = [0 for i in range(self.n)]
self.cplexSolver.define_variable(self.int_var_name_types[i], "integer", objective_coef_vector)
# add constraints
#1. const_names_types
for i in range(self.m):
the_rows = [[self.int_var_name_types[i], [1 for i in range(self.n)]]]
self.cplexSolver.define_constraint([self.const_names_types[i]], the_rows, [self.const_sense_types[i]],
[self.const_rhs_types[i]])
#2 const_name_pattern
coef_list = self.small_rod_size + [-self.large_rod_size[0]]
for j in range(self.n):
temp_var_name_list = []
for i in range(self.m):
temp_var_name_list.append(self.int_var_name_types[i][j])
the_rows = [temp_var_name_list + [self.bin_var_names[j]] ,coef_list]
self.cplexSolver.define_constraint([self.const_name_pattern[j]], [the_rows],
[self.const_sense_pattern[j]], [self.const_rhs_pattern[j]])
def solve_the_model(self):
# Now time to solve it
self.cplexSolver.solve()
def write_the_model(self):
self.cplexSolver.write("cutting_stock_integer.LP")
def get_objective_value(self):
# Get the objective function. - which is number of large-sized rod
return self.cplexSolver.get_objective_value()
if __name__ == "__main__":
if len(sys.argv) != 2:
print("Data is not provided accuratly, solving a cutting stock problem from the default-example")
cs = cutting_stock(DATA.n, DATA.m, DATA.small_rod_demands,DATA.small_rod_size, DATA.large_rod_size)
# Create the model
cs.create_mathematical_model()
if DEBUG:
cs.write_the_model()
# Solve it
cs.solve_the_model()
# obtain the value
print("optimal rod used", int(cs.get_objective_value()))
else:
extract_data_from_input_file = sys.argv[1]
print(extract_data_from_input_file)
```

The solution obtained from the above script is 95 large-sized rods.

The above formulated model was already studied in 1939 by a famous mathematician “L. V. KANTOROVICH”. The paper can be downloaded from the following link: It studied for the first time “Column Generation for Cuttings Stock”

http://resistir.info/livros/kantorovich_mathematical_methods.pdf

Some important points that one can find out from the above model are:

- It is an integer linear program: Difficult problem (NP-hard problem)
- It required many constraints and variables just for the small-sized problem. In our example, for obtaining three different sizes of the small-sized rods, the number of variables generated = n + m*
*n = 100 + 100**3 = 400 ( for the estimate of n, which is 100). - Approximation might not be helpful: Relaxing integer variables to Linear gives us 88 large-sized rods – A lower bound – the gap is 7. For a bigger problem, this gap will be more large.
- A practical limitation: symmetry: obtaining small-sized rods from two different large-sized rods leads to different solutions regarding variables’ values, but absolutely the same in practice. For example, in the current case,
, there will be many alternate solutions to the problem with the same optimal value of 95.

## Can we think of some other way to model it ?

Idea: In the previous model, more than one binaries could be combined if they belong to those large-sized rods that produces the same number/types of small-sized rods. The idea of `Pattern’ hits into our mind, right ?

Now, imagine if we have large-sized rods with various patterns instead of the number of large-sized rods. In that case. The objective is to minimize the number of such different patterns.

Objective function

Here

We also consider a decision variable

Constraints will be:

1.

2.

3.

4.

In the above formulation, if the pattern information is available, the problem become an MIP and can be easily solved. Now, there are two very strong features with this model that make it widely acceptable and benefit us in solving the cutting stock problem efficiently. 1) non-symmetry, i.e., unlike the previous model, a KANTOROVICH model, we do not have different solution set set with same objective value. 2) very strong relaxation: The gap between optimum value of the model and its optimum relaxation is always less than 2 . In fact, there is no instance known with a gap greater than

Now let us discuss the technique, a **column generation technique** which solves this pattern-based formulation-

## Column Generation for Cutting Stock Problem

From the previous paragraph, we got to know a mathematical formulation of the cutting stock problem which is pattern based. But the problem with this formulation is obtaining the patterns as there will be many variables/columns in the model. What, if **we start with few columns and generate the new by some intelligence, if only requires** (the one which gives rise the better objective value to the model)?

The algorithm, **column generation or delayed column generation** does this work. We solve the problem in iterations. In the initial iteration phase, we start with some intuitive valid patterns and solve the problem we call the problem master problem. To compute a better valid pattern we solve a separate problem (subproblem). We also call it a **Knapsack problem**. We will discuss Knapsack problem and its method to solve in other Tutorial.

To understand Column generation well, let us come back to our previous example.

We start with initial iteration:

Iteration 1:

Master:

minimize

Subject to,

To solve, we need to have a valid pattern:[80 50 100] units of [4,6,7] sized rods can be accommodated by [1 0 0], [0 ,1, 0] and [0, 0 , 1] patterns. So in that case. We redefine the model. Master becomes

**Master problem**:

minimize

subject to

All the variables are continuous (assume the problem is linear).

Let us solve them

But still it has many variables, let us drop the variables

We call the resulting model as a **restricted master problem**:

Restricted Master:

minimize

subject to

All the variables are continuous.

The optimal solution to the model is 230. the dual value to the corresponding constraints are =

Now to check the solution is also optimal to the **master problem**. We need to check whether any reduced cost corresponds to the new non basic column

For example, For column

reduced cost,

If we are able to find the value of

Certainly, I can easily find at least one pattern that respects the pattern size restriction.

It is **[2 0 1]**. If we check the reduced cost =

Since, the reduced cost is negative, we have direction of improvement (Are you familiar with this term ? We will discuss in our Tutorial on Simplex algorithm, Later). We will add this pattern/column to the restricted master problem. When solving the updated restricted master problem, the column

For other non-basic columns, we compute the reduced cost,

We select the column (pattern) to be added to the master problem if it leads to negative reduced costs.

Do determine it, we can write a mathematical model

minimize

Subject to,

Where,

If we look this sub-problem carefully, its a special problems of filling up knapsack or a bag of size (say weight) 15 unit by smaller items of different weights. The objective function coefficient associated to the associated variables are the valuation (importance) of each items. We can solve them by mathematical programming solvers or using other heuristic/dynamic programming methods, which mostly far more faster than the IP solver for large number of variables.

Let us implement column generation method and check the optimal cut of rods in our example.

## A basic Column Generation for Cuttings Stock: A python Implementation for our problem

```
# -*- coding: utf-8 -*-
import cplex
import sys
import numpy as np
# Input data
DEBUG = True
class DATA:
n= 100 # number of large-sized rod
m = 3 # number of required rod types
small_rod_demands = [80,50,100]
small_rod_size = [4.,6.,7.]
large_rod_size = [15.0]
class CPLEX:
def __init__(self):
self.model = cplex.Cplex()
self.sub_created = False
def solve_model(self, master_flag):
if master_flag:
self.model.set_problem_type(self.model.problem_type.LP) # you need to specify the type of model, else it will take MIP
else:
self.model.set_problem_type(self.model.problem_type.MILP)
self.model.solve()
def write_problem(self, file_name):
self.model.write(file_name)
def get_objective_value(self):
return self.model.solution.get_objective_value()
def get_solution(self):
return self.model.solution.get_values()
def define_variable(self, varname_list, type_of_var, objective_coef_vector):
n_var = len(varname_list)
if type_of_var == "binary":
self.model.variables.add(obj=objective_coef_vector , names =varname_list,
types = [self.model.variables.type.binary] * n_var )
elif type_of_var == "continuous":
self.model.variables.add(obj= objective_coef_vector, names =varname_list,
types = [self.model.variables.type.continuous] * n_var)
elif type_of_var == "integer":
self.model.variables.add(obj= objective_coef_vector, names =varname_list,
types = [self.model.variables.type.integer] * n_var)
def define_constraint_by_row(self, the_const_names, the_rows, the_senses, my_rhs ):
self.model.linear_constraints.add(lin_expr=the_rows, senses= the_senses, rhs = my_rhs, names = the_const_names)
def get_dual(self):
return self.model.solution.get_dual_values()
def get_pattern(self,pattern_lot,n_const):
return [[self.model.linear_constraints.get_coefficients(i,j) for i in range(n_const)]
for j in range(len(pattern_lot)) if pattern_lot[j]!=0]
def create_sub_generate_pattern(self, coef_multiplier, n_var, large_rod_size, small_rod_size):
if not self.sub_created:
self.var_name = ["y_"+str(i) for i in range(n_var)]
self.model.variables.add(obj=coef_multiplier,types=[self.model.variables.type.integer]*n_var,
names=self.var_name)
row = [[self.var_name,small_rod_size]]
the_senses = 'L'
const_name = 'P'
self.model.linear_constraints.add(lin_expr=row,senses=[the_senses], rhs = large_rod_size,
names= [const_name])
self.model.objective.set_sense(self.model.objective.sense.maximize)
self.sub_created = True
self.model.objective.set_linear([(i,coef_multiplier[i]) for i in range(n_var)])
self.model.solve(False) # For solving a subproblem, Pass False as argument, True is for Master
if self.get_objective_value() > 1: # Since min reduced cost = 1 - dual*variable <0: max dual*variable - 1 >0 => max dual*variable >1
return self.model.solution.get_values()
else:
return []
def update_master(self, pattern,iteration, n_const):
# add a variable to master
self.model.variables.add(obj=[1],types=[self.model.variables.type.continuous],
names=['y_'+str(iteration)], columns=[[[i for i in range(n_const)],pattern]])
class cutting_stock_column_generation:
def __init__(self, n, m, small_rod_demands,small_rod_size,large_rod_size):
self.n = n
self.m = m
self.small_rod_demands = small_rod_demands
self.small_rod_size = small_rod_size
self.large_rod_size = large_rod_size
self.master_flag = True
self.initial_pattern = np.identity(m) # different demand types are m
# initialize restricted master
self.int_var_pattern_count = ["y_"+str(i) for i in range(m)]
self.const_names_types = ["t_"+str(j) for j in range(m)]
self.const_sense_types = ["G" for j in range(m)]
self.const_rhs_types = small_rod_demands
self.cplexSolver1 = CPLEX() # For master
self.cplexSolver2 = CPLEX() # For sub
def create_restricted_master_model(self):
# 1. number of types of patterns
objective_coef_vector = [1 for i in range(self.m)]
self.cplexSolver1.define_variable(self.int_var_pattern_count, "continuous", objective_coef_vector)
# add constraints
#1. const_names_types
for i in range(self.m):
the_rows = [[self.int_var_pattern_count, list(self.initial_pattern[i])]]
self.cplexSolver1.define_constraint_by_row([self.const_names_types[i]],
the_rows, [self.const_sense_types[i]],
[self.const_rhs_types[i]])
def iterative_loop(self, WRITE_DEBUG):
# get the dual values of the initial restricted master:
# solve the initial master
objvalue =0.
dual_value = self.cplexSolver1.get_dual()
pattern = self.cplexSolver2.create_sub_generate_pattern(dual_value, self.m, self.large_rod_size,
self.small_rod_size )
if WRITE_DEBUG:
self.cplexSolver1.write_problem("initial_master.LP")
self.cplexSolver2.write_problem("initial_sub.LP")
column = self.m
while(pattern):
if pattern:
self.cplexSolver1.update_master(pattern,column, self.m)
self.cplexSolver1.solve_model(self.master_flag)
dual_value = self.cplexSolver1.get_dual()
pattern = self.cplexSolver2.create_sub_generate_pattern(dual_value,self.m, self.large_rod_size,
self.small_rod_size )
column = column + 1
if WRITE_DEBUG:
self.cplexSolver1.write_problem("initial_master_"+str(column - self.m)+".LP")
self.cplexSolver2.write_problem("initial_sub_"+str(column - self.m)+".LP")
else:
objvalue = self.cplexSolver1.get_objective_value()
break
return objvalue, pattern
def solve_the_model_master(self):
# Now time to solve it
self.cplexSolver1.solve_model(self.master_flag)
def get_solution(self):
return self.cplexSolver1.get_solution()
def get_objective_value(self):
# Get the objective function. - which is number of large-sized rod
return self.cplexSolver1.get_objective_value()
def get_pattern(self,pattern_lot):
# Get columns
return self.cplexSolver1.get_pattern(pattern_lot, self.m)
if __name__ == "__main__":
WRITE_DEBUG = True
if len(sys.argv) != 2:
print("Data is not provided accuratly, solving a cutting stock problem from the default-example")
cs = cutting_stock_column_generation(DATA.n, DATA.m, DATA.small_rod_demands,DATA.small_rod_size, DATA.large_rod_size)
# Create the model
cs.create_restricted_master_model()
#if DEBUG:
# cs.write_the_model_master("initial_restricted_master.LP")
cs.solve_the_model_master()
cs.iterative_loop(WRITE_DEBUG)
# obtain the value
print("optimal rod used", int(cs.get_objective_value()))
pattern_lot = cs.get_solution()
print("optimal pattern lot", cs.get_solution())
print("optimal pattern", cs.get_pattern(pattern_lot))
else:
extract_data_from_input_file = sys.argv[1]
print(extract_data_from_input_file)
```

The optimal solution we get is 95 rods and the patterns are: [[2.0, 0.0, 1.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]]

## Fork the above code from my git repository for your use.

## Other benefit in using Column generation technique and its application apart from Cutting stock problem

- Obtain lower and upper bounds of the pattern of large-sized rod used and so its partial solution without completely solving the problem
- Adapting heuristic in the solving process – knapsack problem can use any dynamic/heuristic method for speeding up the solving time
- For most of the mathematical models where number of columns are very large, column-generation can be the preferred option to try.