Description: This notebook describes how to implement column generation, which is a large scale optimization scheme, in JuMP. The cutting stock problem has been used as an illustrative example.

Author: Shuvomoy Das Gupta

License: Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Using Julia+JuMP for optimization - column generation


Implementing large scale optimization techniques such as column generation is really easy using JuMP. To explain how to implement column generation in JuMP, we consider the famous cutting stock problem. For more details about the problem, see pages 234-236 of Introduction to Linear Optimization by Bertsimas and Tsitsiklis.

Notation and notions:

  • Width of a large roll is $W$, and it needs to be cut into smaller width papers according to customer demand
  • The set of indices of all feasible patterns is, $\mathcal{J}=\{1,2,\ldots,n\}$, where $n$ is a very large number
  • A strict subset of $\mathcal{J}$ that is considered in the master problem is $\mathcal{J}'$
  • The dummy index for a pattern is $j$
  • The index set of all possible paper-widths is, $\mathcal{M}=\{1,2,\ldots,m\}$
  • The width of the paper with index $i$ is $w_i$
  • The demand for the paper of width $w_i$ is $b_i$
  • Number of smaller rolls of width $w_i$ produced by pattern $j$ is denoted by $a_{ij}$
  • Number of large rolls cut according to pattern $j$ is denoted by $x_j$

Original unabridged problem:

$$ \begin{align} &\text{minimize} && \sum_{j \in \mathcal{J}}{x_j} \\ &\text{subject to} &&\\ & &&\forall i \in \mathcal{M} \quad \sum_{j \in \mathcal{J}}{a_{ij} x_j}=b_i \\ & && \forall j \in \mathcal{J} \quad x_j \geq 0 \\ \end{align} $$

Because the set $\mathcal{J}$ can be astronomically large, even storing the problem is a challenge. So, we start with a smaller version of the problem, called the master problem, by replacing $\mathcal{J}$ with a strict subset $\mathcal{J}'$, which is much smaller than the original one.


Structure of the decomposition

Master Problem: $$ \begin{align} &\text{minimize} && \sum_{j \in \mathcal{J}'}{x_j} \\ &\text{subject to} &&\\ & &&\forall i \in \mathcal{M} \quad \sum_{j \in \mathcal{J}'}{a_{ij} x_j}=b_i \\ & && \forall j \in \mathcal{J}' \quad x_j \geq 0 \\ \end{align} $$

After solving the master problem, we want to check the optimality status. Structure of the cutting stock problem allows us to construct a subproblem which can do this very easily.

Subproblem: $$ \begin{align} &\text{minimize} && 1 - \sum_{i \in \mathcal{M}} \quad {p_i a_{i {j^*}}} \\ &\text{subject to} &&\\ & && \forall i \in \mathcal{M} \quad a_{i {j^*}} \geq 0, \quad a_{ij^*} \; \text{integer} \\ & && \sum_{i \in \mathcal{M}}{w_i a_{i{j^*}}} \leq W\\ \end{align} $$

The objective of the subproblem is the minimum of the reduced cost vector of the original problem. If the objective value of the subproblem is greater than or equal to $0$, then the current solution of the master problem is optimal for the original unabridged problem. Otherwise, add the resultant cost reducing column $(a_{i {j^*}})_{i \in \mathcal{M}}=A_{j*}$ and a corresponding new variable $x_{j*}$ is added to the master problem. The modified master problem is as follows:

Modified Master Problem $$ \begin{align} &\text{minimize} && \sum_{j \in \mathcal{J}'}{x_j} + x_{j^*} \\ &\text{subject to} &&\\ & &&\forall i \in \mathcal{M} \quad \sum_{j \in \mathcal{J}'}{a_{ij} x_j}+a_{i j^*} x_{j^*}=b_i \\ & && \forall j \in \mathcal{J}' \quad x_j \geq 0, x_j^* \geq 0 \\ \end{align} $$

The pseudocode for the cutting stock problem is given below.

Pseduocode

  • Input preliminary data for starting the problem
  • Solve the master problem with the initial data
$$ \begin{align} &\text{minimize} && \sum_{j \in \mathcal{J}'}{x_j} \\ &\text{subject to} &&\\ & &&\forall i \in \mathcal{M} \quad \sum_{j \in \mathcal{J}'}{a_{ij} x_j}=b_i \\ & && \forall j \in \mathcal{J}' \quad x_j \geq 0 \\ \end{align} $$
  • Collect the dual variables for the equality constraints and store them in an array $(p_i)_{i \in \mathcal{M}}$

  • Solve the sub problem

$$ \begin{align} &\text{minimize} && 1 - \sum_{i \in \mathcal{M}} \quad {p_i a_{i {j^*}}} \\ &\text{subject to} &&\\ & && \forall i \in \mathcal{M} \quad a_{i {j^*}} \geq 0, \quad a_{ij^*} \; \text{integer} \\ & && \sum_{i \in \mathcal{M}}{w_i a_{i{j^*}}} \leq W\\ \end{align} $$
  • Flow control:

while ( $\text{optimal value of the subproblem} < 0$)

  • Add the column $(a_{i {j^*}})_{i \in \mathcal{M}}=A_{j*}$ to $A$
  • Add a corresponding new variable $x_{j*}$ to the list of variables
  • Solve the modified master problem

$$ \begin{align} &\text{minimize} && \sum_{j \in \mathcal{J}'}{x_j} + x_{j^*} \\ &\text{subject to} &&\\ & &&\forall i \in \mathcal{M} \quad \sum_{j \in \mathcal{J}'}{a_{ij} x_j}+a_{i j^*} x_{j^*}=b_i \\ & && \forall j \in \mathcal{J}' \quad x_j \geq 0 \\ & && \qquad \qquad \; \; x_{j^*} \geq 0 \end{align} $$

  • Collect the dual variables for the equality constraints and store them in an array $(p_i)_{i \in \mathcal{M}}$
  • Solve the sub problem as before
  • Set $\mathcal{J}':=\mathcal{J}'\cup \{j^*\}$

end while

  • Display the results

Master Problem Modification in JuMP

The problem modification can be done by using the already mentioned @variable macro:

$$ \texttt{@variable}(m, l \leq x_\text{new} \leq u, \texttt{Int}, \texttt{objective} = c_\text{new}, \texttt{inconstraints} = \text{arrayConstrrefs}, \texttt{coefficients} = \text{arrayCoefficients}) $$

Here:

  • The name of the original model is $m$.
  • The new variable to be added is $x_\text{new}$ with lower bound $l$ and upper bound $u$.
  • The type of the variable can be Int, Bin. For real variable the third argument is left vacant.
  • The original objective, say $f_o(x)$ will become $f_o(x) + c_\text{new} x_\text{new}$ after modification
  • The array $\texttt{arrayConstrrefs}$ contain references to those constraints that need to be modified by inclusion of $x_\text{new}$
  • The array $\texttt{arrayCoefficients}$ contain the coefficients that have to multiplied with $x_\text{new}$ and then added to the constraints referenced by $\texttt{arrayConstrrefs}$ in an orderly manner. For example, if the $i$th element of $\texttt{arrayConstrrefs}$ refers to a constraint $a_i^T x \lesseqgtr b_i$, then after invoking the command, the constraint is modified as: $a_i^T x +\texttt{arrayCoefficients}[i] x_\text{new} \lesseqgtr b_i$

Implementing one iteration of the column generation algorithm

To understand how the column generation is working in Julia, we implement one iteration of the column generation algorithm manually. The entire code is presented in the next section.

In [1]:
# Uploading the packages:
# -----------------------

using JuMP 

# We will use default solvers
In [2]:
# Input preliminary data for starting the problem
# -----------------------------------------------

W=100
cardinalityM=5
M=collect(1:cardinalityM)
A=eye(cardinalityM)
p=zeros(5)
b=[45; 38; 25; 11; 12]
w=[22; 42; 52; 53; 78]
Out[2]:
5-element Array{Int64,1}:
 22
 42
 52
 53
 78
In [3]:
# Description of the master problem with the initial data
#----------------------

cutstockMain = Model() # Model for the master problem
Jprime=collect(1:size(A,2)) # Initial number of variables
@variable(cutstockMain, 0 <= x[Jprime] <= 1000000) # Defining the variables
@objective(cutstockMain, Min, sum(1*x[j] for j in Jprime)) # Setting the objective
@constraint(cutstockMain, consRef[i=1:cardinalityM], sum(A[i,j]*x[j] for j in Jprime)==b[i]) 
# Here the second argument consRef[i=1:cardinalityM] means that the i-th constraint aᵢᵀx = bᵢ has the corresponding constraint reference
# consRef[i]
print(cutstockMain)
Min x[1] + x[2] + x[3] + x[4] + x[5]
Subject to
 x[1] == 45
 x[2] == 38
 x[3] == 25
 x[4] == 11
 x[5] == 12
 0 <= x[i] <= 1.0e6 for all i in {1,2,..,4,5}
In [4]:
# Solving the master problem with the initial data
# ------------------------------------------------
solve(cutstockMain)
println("Current solution of the master problem is ", getvalue(x))
println("Current objective value of the master problem is ", getobjectivevalue(cutstockMain))
Current solution of the master problem is x: 1 dimensions:
[1] = 45.0
[2] = 38.0
[3] = 25.0
[4] = 11.0
[5] = 12.0

Current objective value of the master problem is 131.0
In [5]:
#Collect the dual variables for the equality constraints and store them in an array p
for i in M
    p[i] = getdual(consRef[i]) # These p[i] are the input data for the subproblem
end 
println("The array storing the dual variables is ", p)
The array storing the dual variables is [1.0,1.0,1.0,1.0,1.0]
In [6]:
# Describe the sub problem
# ------------------------
cutstockSub=Model() # Model for the subproblem
@variable(cutstockSub, 0 <= Ajstar[M] <= 1000000, Int )
@objective(cutstockSub, Min, 1-sum(p[i]*Ajstar[i] for i in M))
@constraint(cutstockSub, sum(w[i]*Ajstar[i] for i in M) <= W)
print(cutstockSub)
Min -Ajstar[1] - Ajstar[2] - Ajstar[3] - Ajstar[4] - Ajstar[5] + 1
Subject to
 22 Ajstar[1] + 42 Ajstar[2] + 52 Ajstar[3] + 53 Ajstar[4] + 78 Ajstar[5] <= 100
 0 <= Ajstar[i] <= 1.0e6, integer, for all i in {1,2,..,4,5}
In [7]:
# Solve the sub problem
# ---------------------
solve(cutstockSub)
minreducedCost=getobjectivevalue(cutstockSub)
println("The minimum component of the reduced cost vector is ", minreducedCost)
The minimum component of the reduced cost vector is -3.0

The minimum component of the reduced cost vector is negative, so we have a suboptimal solution.

In [8]:
if minreducedCost >= 0
    println("We are done, current solution of the master problem is optimal")
else
    println("We have a cost reducing column ", getvalue(Ajstar))
end
We have a cost reducing column Ajstar: 1 dimensions:
[1] = 4.0
[2] = 0.0
[3] = 0.0
[4] = 0.0
[5] = 0.0

In [9]:
typeof(Ajstar)
Out[9]:
JuMP.JuMPArray{JuMP.Variable,1,Tuple{Array{Int64,1}}}

Now Ajstar is of type JuMPDict. To use it in the modified master problem, we have to store values from Ajstar in a column vector.

In [10]:
Anew=Float64[] # This Anew correspond to the newly added column to the A matrix
for i in 1:cardinalityM
    push!(Anew, getvalue(Ajstar)[i])
end

When we add the cost reducing column Anew to the original matrix A, it also gives rise to a new variable xNew corresponding to Anew. Now we want to keep track of the new variables that are added by the subproblem. We do this by declaring an array of Variables named xNewArray, which will contain all such newly added variables in the process of column generation.

In [11]:
xNewArray=Variable[] # The newly added variables by flow control will be
# pushed to the new array of variables xNewArray
Out[11]:
$$ Empty Array{Variable} (no indices) $$

Here we just illustrate one iteration of the while loop manually, because, for now, we are interested to understand how JuMP is managing the flow control and modifying the master problem and the sub problem.

Let's modify the master problem by adding the new column Anew to the old A matrix. Note that we do not have to rewrite the entire model.

In [12]:
# Modify the master problem by adding the new column Anew to the old A matrix
@variable(
cutstockMain, # Model to be modified
0 <= xNew <= 1000000, # New variable to be added
objective=1, # cost coefficient of new variable in the objective
inconstraints=consRef,  # constraints to be modified
coefficients=Anew # the coefficients of the variable in those constraints
) 

# The line above adds the column (aᵢⱼ*)ᵢ=Aⱼ* to A <br>
# and add a corresponding new variable xⱼ* to the list of variable

push!(xNewArray, xNew) # Pushing the new variable in the array of new variables
print(cutstockMain)
Min x[1] + x[2] + x[3] + x[4] + x[5] + xNew
Subject to
 x[1] + 4 xNew == 45
 x[2] == 38
 x[3] == 25
 x[4] == 11
 x[5] == 12
 0 <= x[i] <= 1.0e6 for all i in {1,2,..,4,5}
 0 <= xNew <= 1.0e6

Though we are showing only one iteration of the flow control, in the final code for sure we want to have a

while ( some condition )
(
...
)
end

block.

Now if we do not do anything else in the final code, all the names of the newly added variables by the while loop will be the same: xNew! JuMP is intelligent enough to treat them as separate variables, but it is not very human-friendly. It is more convenient if the newly added variables were given different names, which we can achieve by setName(oldName, newName) function.

In [13]:
setname(xNew, string("x[",size(A,2)+1,"]")) # Changing the name of the variable 
# otherwise all the newly added variables will have name xNew!
# size(A,2) gives the column number of A
Out[13]:
"x[6]"

Let us see if the name of the variable has changed as desired.

In [14]:
print(cutstockMain) # Let us see if the name of the variables have changed as desired
Min x[1] + x[2] + x[3] + x[4] + x[5] + x[6]
Subject to
 x[1] + 4 x[6] == 45
 x[2] == 38
 x[3] == 25
 x[4] == 11
 x[5] == 12
 0 <= x[i] <= 1.0e6 for all i in {1,2,..,4,5}
 0 <= x[6] <= 1.0e6

Indeed it has! Now let's solve the modified master problem, and then collect the associated dual variables for the equality constraints and store them in the array p.

In [15]:
statusControlFlow=solve(cutstockMain) # Solve the modified master problem

getdual(consRef)
for i in M
    p[i] = getdual(consRef)[i] 
end 

println(p)
[0.25,1.0,1.0,1.0,1.0]

Now we solve the subproblem for the current solution of the master problem:

In [16]:
# Solving the modified sub problem        
@variable(cutstockSub, 0 <= Ajstar[M] <= 1000000, Int )
@objective(cutstockSub, Min, 1-sum(p[i]*Ajstar[i] for i in M))
@constraint(cutstockSub, sum(w[i]*Ajstar[i] for i in M) <= W)
print(cutstockSub) # Let's see what is the current subproblem looks like
solve(cutstockSub)
minReducedCost=getobjectivevalue(cutstockSub)
println("Current value of the minimum of the reduced cost vector is ", minReducedCost)
Min -0.25 Ajstar[1] - Ajstar[2] - Ajstar[3] - Ajstar[4] - Ajstar[5] + 1
Subject to
 22 Ajstar[1] + 42 Ajstar[2] + 52 Ajstar[3] + 53 Ajstar[4] + 78 Ajstar[5] <= 100
 22 Ajstar[1] + 42 Ajstar[2] + 52 Ajstar[3] + 53 Ajstar[4] + 78 Ajstar[5] <= 100
 0 <= Ajstar[i] <= 1.0e6, integer, for all i in {1,2,..,4,5}
 0 <= Ajstar[i] <= 1.0e6, integer, for all i in {1,2,..,4,5}
Current value of the minimum of the reduced cost vector is -1.0
WARNING: Solver does not appear to support providing initial feasible solutions.

The optimal value of the current subproblem is negative (which will be tested by the conditional statement of the while loop in the final code), giving us a cost reducing column to be added in the master problem. As before we have to store the column Ajstar in a column vector Anew.

In [17]:
#Store the components of the solution of current subproblem into the column Anew    
Anew=Float64[]
for i in 1:cardinalityM
    push!(Anew, getvalue(Ajstar)[i])
end

println("New column to be added to A is: ", Anew)
New column to be added to A is: [0.0,2.0,0.0,0.0,0.0]

Okay, we have understood how JuMP is working in the column generation process. The entire code of the cutting stock problem is given below:

Cutting stock problem code:

In [18]:
# Verfied to be working:

# Uploading the packages:
# -----------------------

using JuMP
using GLPKMathProgInterface

# Input preliminary data for starting the problem
# -----------------------------------------------

W=100
cardinalityM=5
M=collect(1:cardinalityM)
A=eye(cardinalityM)
p=zeros(5)
b=[45; 38; 25; 11; 12]
w=[22; 42; 52; 53; 78]

@time begin # time measurement begins

# Solve the master problem with the initial data
#-----------------------------------------------

cutstockMain = Model() # Model for the master problem
Jprime=collect(1:size(A,2)) # Intial number of variables
@variable(cutstockMain, 0 <= x[Jprime] <= 1000000) # Defining the variables
    
@objective(cutstockMain, Min, sum(1*x[j] for j in Jprime)) # Setting the objective
    
@constraint(cutstockMain, consRef[i=1:cardinalityM], sum(A[i,j]*x[j] for j in Jprime)==b[i]) # Adding the constraints
# Here the second argument consRef[i=1:cardinalityM] means that the i-th constraint aᵢᵀx = bᵢ has 
# the corresponding constraint reference consRef[i]

solve(cutstockMain)
 
#Collect the dual variables for the equality constraints and store them in an array p
getdual(consRef)
for i in M
    p[i] = getdual(consRef)[i] # These p[i] are the input data for the subproblem
end 

# Solve the sub problem
# -------------------

cutstockSub=Model() # Model for the subproblem
@variable(cutstockSub, 0 <= Ajstar[M] <= 1000000, Int )
@objective(cutstockSub, Min, 1-sum(p[i]*Ajstar[i] for i in M))
@constraint(cutstockSub, sum(w[i]*Ajstar[i] for i in M) <= W)
solve(cutstockSub)
minReducedCost=getobjectivevalue(cutstockSub)

Anew=Float64[] # This Anew correspond to the newly added column to the A matrix
for i in 1:cardinalityM
    push!(Anew, getvalue(Ajstar)[i])
end

xNewArray=Variable[] # The newly added variables by flow control will be pushed to the new array of variables xNewArray

k=1 # Counter for the while loop

                                                # Flow control
                                                # ------------

while minReducedCost < 0 #while (current solution of the master problem is suboptimal, i.e., subproblem objective value < 0)
    # Solve the master problem by adding the new column Anew to the old A matrix
     @variable(
             cutstockMain, # Model to be modified
             0 <= xNew <= 1000000, # New variable to be added
             objective=1, # cost coefficient of new varaible in the objective
             inconstraints=consRef,  # constraints to be modified
             coefficients=Anew # the coefficients of the variable in those constraints
            ) 
        # The line above adds the column (aᵢⱼ*)ᵢ=Aⱼ* to A <br>
        # and add a corresponding new variable xⱼ* to the list of variable
        push!(xNewArray, xNew) # Pushing the new variable in the array of new variables
        setname(xNew, string("x[",size(A,2)+k,"]")) # Changing the name of the variable 
                                               # otherwise all the newly added variables will have name xNew!
        k=k+1 # Increasing k by 1
        statusControlFlow=solve(cutstockMain)

        #Collect the dual variables for the equality constraints and store them in an array p
        getdual(consRef)
        for i in M
            p[i] = getdual(consRef)[i] 
        end 

        # Solving the modified sub problem        
        @variable(cutstockSub, 0 <= Ajstar[M] <= 1000000, Int )
        @objective(cutstockSub, Min, 1-sum{p[i]*Ajstar[i],i in M})
        @constraint(cutstockSub, sum{w[i]*Ajstar[i], i in M} <= W)
        solve(cutstockSub)
        minReducedCost=getobjectivevalue(cutstockSub)

        #Store the components of the solution of current subproblem into the column Anew    
        Anew=Float64[]
        for i in 1:cardinalityM
            push!(Anew, getvalue(Ajstar)[i])
        end
    end # While loop ends
    
end # time measurement ends

# Print the results
# -----------------

println("Objective value: ", getobjectivevalue(cutstockMain))
println("Current Solution is: ", getvalue(x))
println("With ", length(xNewArray), " variables added by flow control:")
for i in 1:length(xNewArray)
    println("[",size(A,2)+i,"] = ",getvalue(xNewArray[i]))
end
println("Reduced cost of the current solution is ", getobjectivevalue(cutstockSub))
  0.003651 seconds (4.00 k allocations: 265.625 KB)
Objective value: 57.25
Current Solution is: x: 1 dimensions:
[1] = 0.0
[2] = 0.0
[3] = 0.0
[4] = 0.0
[5] = 0.0

With 5 variables added by flow control:
[6] = 8.25
[7] = 1.0
[8] = 11.0
[9] = 25.0
[10] = 12.0
Reduced cost of the current solution is 0.0