Introduction to large scale programming: column generation
¶

Alfred Galichon (NYU & Sciences Po) and Antoine Jacquet (Sciences Po)
¶

'math+econ+code' masterclass series
¶

With python code examples
¶

© 2018–2026 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed here.

If you reuse material from this masterclass, please cite as:
Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/

Learning objectives¶

  • Large-scale linear programming
  • Column generation
  • The cutting stock problem
  • The knapsack problem

References¶

  • Dantzig & Wolfe (1960). Decomposition Principle for Linear Programs. Operations Research

Libraries¶

Motivation¶

In this lecture, we introduce large-scale linear programming through the method of column generation, using the cutting stock problem as our main motivating example and the knapsack problem as its key computational building block. In the cutting stock problem, each feasible cutting pattern can be viewed as a column in a very large linear program, typically far too large to enumerate in full. Column generation addresses this difficulty by solving a restricted master problem and then adding new columns only when they are useful. The crucial step is the pricing problem: given the current dual values, finding a new improving cutting pattern turns out to be a knapsack problem. Studying these two problems together therefore provides a natural and concrete introduction to both the modeling logic of large-scale LP and the algorithmic structure of column generation.

The cutting stock problem¶

A paper maker produces raw rolls of length $W \in \mathbb R_+$.

The rolls are cut into final products of various types $i\in[I]$, which are characterized by their length $w_i \in \mathbb R_+$. The demand for product $i$ is $q_i$.

A configuration $a\in \mathcal A \subseteq \mathbb N^I$ is a way to cut a raw roll into final products, where $a_{i}$ is the number of times product $i$ appears in configuration $a\in \mathcal A$. As the sum of the lengths of the final products needs to be less or equal than the length of a roll, we need to have
$\sum_{i\in [I]} a_{i} w_i \leq W,$
and, as a result, the cardinality of the set of feasible configurations $\mathcal A = \{ a \in \mathbb N ^I: \sum_{i\in [I]} a_{i} w_i \leq W \}$ can be a huge number.

We assume that each roll has a cost of one no matter in which configuration it is cut, and we let $x(a)$ be the number of rolls cut in a configuration $a \in \mathcal A$. The problem is to determine the number of rolls needed in various configurations to meet the demand at minimal cost. This is an integer programming problem, but we shall drop the requirement that $x(a)$ should be integral and study its linear programming relaxation:

$\min_{x \in \mathbb R^{\mathcal A}_+} \sum_{a\in \mathcal A} x(a)$
s.t. $\sum_{a\in \mathcal A} a_{i} x(a) = q_i,~\forall i \in [I].$

This problem is a difficult problem because the number of variables $|\mathcal A|$ is potentially huge.

Let's take an example:

In [1]:
import numpy as np
w_i = np.array([1365, 1545, 1510, 1735, 1805, 1895, 1915, 2015, 2040, 2125, 2110, 2175, 2245])
q_i = np.array([23, 27, 11, 16, 19, 17, 22, 9, 14, 13, 18, 20, 19])
W = 4500
I = len(w_i)

Let us give a first example of a configuration $(a^{ex}_i)$ with one order $i=0$, one $i=1$, and one $i=2$, and let's test its feasibility:

In [2]:
aex_i = np.array([1,1,1]+[0]*(I-3))
print('aex_i',aex_i)
print("Feasible" if aex_i @ w_i <= W else "Not Feasible")

aex_i [1 1 1 0 0 0 0 0 0 0 0 0 0]
Feasible

We are going to keep track of a basis $\mathcal B$ of elements of $\mathcal A$ that form a basis of $\mathbb R^I$. $\mathcal B$ has $I$ elements, and we denote by $a^1,...,a^I$ its elements, which are all $I$-dimensional vectors. Let $B=[a^1,..,a^I]$ be the $I\times I$ matrix whose columns are the $a^i$'s, so that the $ij$-th entry of $B$ is $a_i^j$.

In [3]:
B_i_j = np.diag( W // w_i)
print(B_i_j)

[[3 0 0 0 0 0 0 0 0 0 0 0 0]
 [0 2 0 0 0 0 0 0 0 0 0 0 0]
 [0 0 2 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 2 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 2 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 2 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 2 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 2 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 2 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 2 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 2 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 2 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 2]]

Our initial feasible basis will consist in configurations that produce only one type of final product, in the largest possible quantity that one roll may produce, that is for $i \in[I]$, $a^i_i = \lfloor W / w_i \rfloor$ and $a_i^j=0$ for $j \in [I] \backslash \{i\} $.
The basic solution then consists in providing $x_j$ amount of configuration $a^j$. We need to have $\sum_{j \in [I]} a^j_i x_j = q_i$, that is, in matrix terms, $B x = q$,
thus $x = B^{-1} q$, that is, the initial basic solution is
$x^0_j=q_j / \lfloor W / w_j \rfloor$ for $j\in[I]$.

In [4]:
x0_j = np.linalg.solve(B_i_j,q_i)
x0_j
Out[4]:
array([ 7.66666667, 13.5       ,  5.5       ,  8.        ,  9.5       ,
        8.5       , 11.        ,  4.5       ,  7.        ,  6.5       ,
        9.        , 10.        ,  9.5       ])

Assuming in this example that $C(a)=1$ for all $a \in \mathcal A$ (so that $c=0$) the associated objective value is $\sum_{j\in [I]} x^0_j$, that is:

In [5]:
x0_j.sum()
Out[5]:
110.16666666666667

Now let's go back to our example configuration aex_i (which, as we recall, has one each $i=0,1$ and $2$), and let's check if we are making a saving by introducing it in order to replace vectors of the basis.

If we introduce one unit of the example configuration, we need to remove $z_j$ of each of the basic configurations $a^j$, where
$a^{ex}_i = \sum_{j\in [I]} a^j_i z_j$.

This can be expressed as
$a^{ex} = B z,$
so we have
$z = B^{-1} a^{ex}.$

The marginal benefit of this replacement is $\sum_{j\in[I]} z_j - 1 = \mathbf{1}_I^\top B^{-1} a^{ex} -1 = k^\top a^{ex} -1$, where we have let
$k = (B^{-1})^\top \mathbf{1}_I$.

Exercise. Show that the result is not modified if the cost associated with a roll cut in configuration $a\in \mathcal A$ is $1 + c^\top a = 1 + \sum_{i \in [I]} c_i a_i$, where $c\in \mathbb R^I_+$.

In our example, we get:

In [6]:
k_i = np.linalg.solve(B_i_j.T,np.ones(I) )
print("Marginal benefit=",k_i.dot(aex_i )- 1)

Marginal benefit= 0.33333333333333326

This means that introducing this configuration would have a net benefit. According to the standard simplex method, we should let this configuration enter the basis, and determine which configuration exits. In order to do so, we shall introduce a quantity $\epsilon>0$ of the new configuration $a^{ex}$. To make up for it, we need to adjust the quantities of the configurations $a^j$ -- let $x^\epsilon_j$ be the quantity of configuration $a^j$ held. We have $\epsilon a^{ex}+ B x^{\epsilon} = q$, and thus
$x^{\epsilon} = B^{-1} q - \epsilon B^{-1} a^{ex} = x^0 - \epsilon z$. We will set the value of $\epsilon$ as the largest value such that $x^{\epsilon}\geq 0$, that is, we shall take
$\epsilon = \min_{j\in [I]} \{ \frac {x^0_j} {z_j}: z_j>0 \},$ and the corresponding argmin $j^{exit}$ is the index of the column corresponding to the departing configuration.

We then do a standard pivot step and we replace the column

In [7]:
def pivot(B_i_j,aent_i ,xB_j):
    z_j = np.linalg.solve(B_i_j,aent_i)
    thedic = {j: xB_j[j] / z_j[j] for j in range(I) if z_j[j]>0}
    jexit = min(thedic, key = thedic.get)
    epsilon = thedic[jexit]
    xB_j = xB_j - epsilon * z_j
    xB_j[jexit] = epsilon
    B_i_j[:,jexit] = aent_i
    return jexit,xB_j

jexit,x0_j = pivot(B_i_j,aex_i ,x0_j)
print( "Departing column=", jexit, "; new configuration basis matrix=")
print(B_i_j)

Departing column= 2 ; new configuration basis matrix=
[[3 0 1 0 0 0 0 0 0 0 0 0 0]
 [0 2 1 0 0 0 0 0 0 0 0 0 0]
 [0 0 1 0 0 0 0 0 0 0 0 0 0]
 [0 0 0 2 0 0 0 0 0 0 0 0 0]
 [0 0 0 0 2 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 2 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 2 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 2 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 2 0 0 0 0]
 [0 0 0 0 0 0 0 0 0 2 0 0 0]
 [0 0 0 0 0 0 0 0 0 0 2 0 0]
 [0 0 0 0 0 0 0 0 0 0 0 2 0]
 [0 0 0 0 0 0 0 0 0 0 0 0 2]]

At this stage, we wonder which next configuration should enter. When in the first step we chose $a^{ex}$ as an entering configuration, it was arbitrary, for illustration purposes. In practice, we want to find an efficient and systematic way to choose the entering configuration. We look for a configuration that has a positive benefit by computing:
$\max_{a\in \mathcal A} \{ \sum_i k_i a_i \}$
s.t. $\sum_{i\in [I] } a_i w_i \leq W$
$a_i$ integer.
If the value of that problem is $\leq 1$, then no configuration has a positive benefit, and the current base of configuration yields an optimal solution. Otherwise, we shall pick the maximizer as an entering configuration.

This problem is a well-studied problem in combinatorial optimization called the knapsack problem.

The knapsack problem¶

Let $F_I(W)$ be the value of the above problem we make the dependence in $I$ and $W$ explicit. Trying to find the optimal combination of configurations among $i \in [I]$ of total length less than $W$ can be decomposed in finding $a_{I}$, the optimal number configuration of index $I$ and recursively solving the problem among the $i\in[I-1]$ remaining configuration, of total length less than $W-a_{I} w_{I}$. We have the following dynamic programming principle:
$F_I(W) = \max_{a_{I} \in \{ 0,..., \lfloor v / w_{I} \rfloor \} }\{k_{I}a_{I}+F_{I-1}(W-w_{I}a_{I}) \}$
which is initialized by
$F_1(v) = \lfloor v / w_{1} \rfloor \max(k_1,0)$.

In [8]:
def knapsack(k_i,W):
    I=len(k_i)
    F = np.zeros( (I,W+1))
    A = np.zeros( (I,W+1),dtype = int)
    A[0,:]=np.array([(v // w_i[0])  for v in range(W+1)])
    F[0,:]=np.array([(v // w_i[0])* max(k_i[0],0)  for v in range(W+1)])
    for i in range(I-1):
        for v in range(W+1):
            thedic = {a: k_i[i+1]*a +F[i,v-a*w_i[i+1]] for a in range( (v// w_i[i+1])+1 )}
            A[i+1,v] = max(thedic, key = thedic.get)
            F[i+1,v] = thedic[A[i+1,v]]
    v = W
    a_i = np.zeros(I)
    for i in reversed(range(I)):
        a_i[i] = A[i,v]
        v = v-A[i,v]*w_i[i]
    return(a_i)

We compute a candidate entering column $a^{ent}_i$. We check that the candidate entering configuration is indeed feasible, and we check that the candidate entering configuration $a^{ent}$ indeed improves on the cost.

In [9]:
aent_i=knapsack(k_i,W)
print('aent_i=',aent_i)
print("Feasible" if aent_i @ w_i <= W else "Not Feasible")
print("Improving" if 1 < k_i.dot(aent_i) else "Not improving")

aent_i= [1. 2. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Feasible
Improving

If we introduce one unit of the new configuration, we need to remove $z_j$ of each of the basic configurations $j\in B$, where
$a^{ent}_i = \sum_{j\in B} A_{ij} z^{ent}_j$
so we have
$z^{ent} = A^{-1} a^{ent}.$

We look for $\epsilon$ such that
$\sum_{j \in B} A_{ij}(x_j - \epsilon z^{ent}_j)+\epsilon a^{ent}_i$
has one zero term (which defines the $j$ leaving the basis).
This leads us to determine
$\epsilon = min\{ x_j / z_j : z_j >0 \}$.

In [10]:
zent_j = np.linalg.solve(B_i_j,aent_i)
thedic = {j: x0_j[j] / zent_j[j] for j in range(I) if zent_j[j]>0}
jexit = min(thedic, key = thedic.get)
epsilon = thedic[jexit]
epsilon
Out[10]:
8.0

We update the feasible solution and the basis accordingly:

In [11]:
x1_j = x0_j - epsilon * zent_j
x1_j[jexit] = epsilon
print(x1_j)
B_i_j[:,jexit] = aent_i

[ 1.33333333  8.         11.          8.          9.5         8.5
 11.          4.5         7.          6.5         9.         10.
  9.5       ]

We check that the new solution is indeed feasible:

In [12]:
print('xB_j=\n',x1_j)
print('\nslackness=',B_i_j @ x1_j - q_i)

xB_j=
 [ 1.33333333  8.         11.          8.          9.5         8.5
 11.          4.5         7.          6.5         9.         10.
  9.5       ]

slackness= [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

The revised objective value is now:

In [13]:
x1_j.sum()
Out[13]:
103.83333333333334

The column generation algorithm¶

We iterate the previous process and we get:

In [14]:
cont = True
B_i_j = np.diag( W // w_i)
xB_j = np.linalg.solve(B_i_j,q_i)
iter = 0
while cont:
    iter += 1
    print('Obj=',xB_j.sum())
    k_i = np.linalg.solve(B_i_j.T,np.ones(I))
    aent_i=knapsack(k_i,W)
    if 1 >= k_i.dot(aent_i):
        cont=False
    else:
        jexit,xB_j = pivot(B_i_j,aent_i ,xB_j)
        
print( "Converged in "+str(iter)+" steps.\nObjective="+str(xB_j.sum()))
    

Obj= 110.16666666666667
Obj= 105.66666666666667
Obj= 103.83333333333334
Obj= 103.5
Converged in 4 steps.
Objective=103.5