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:
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
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.
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.
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.
Collect the dual variables for the equality constraints and store them in an array $(p_i)_{i \in \mathcal{M}}$
Solve the sub problem
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
The problem modification can be done by using the already mentioned @variable
macro:
Here:
Int
, Bin
. For real variable the third argument is left vacant.$a_i^T x +\texttt{arrayCoefficients}[i] x_\text{new} \lesseqgtr b_i$
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.
# Uploading the packages:
# -----------------------
using JuMP
# We will use default solvers
# 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]
5-element Array{Int64,1}: 22 42 52 53 78
# 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}
# 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
#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]
# 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}
# 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.
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
typeof(Ajstar)
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.
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 Variable
s named xNewArray
, which will contain all such newly added variables in the process of column generation.
xNewArray=Variable[] # The newly added variables by flow control will be
# pushed to the new array of variables xNewArray
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.
# 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.
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
"x[6]"
Let us see if the name of the variable has changed as desired.
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
.
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:
# 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
.
#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:
# 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