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

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

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.

- Input preliminary data for starting the problem
- Solve the master problem with the initial data

Collect the dual variables for the equality constraints and store them in an array $(p_i)_{i \in \mathcal{M}}$

Solve the sub problem

- 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

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

macro:

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$

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

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

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

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

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

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

In [9]:

```
typeof(Ajstar)
```

Out[9]:

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

s 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]:

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

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

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

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

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

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

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

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