【0基础运筹学】【超详细】【附代码】列生成（Column Generation）

列生成（Column generation）是一种解决大型线性程序的有效算法。

相关教程

• 【从零开始】coin-or/CoinUtils Osi Clp Cgl Cbc源码构建debug（CLion/CMake）

相关文献

• Column generation（wikipedia）
• 干货 | 10分钟带你彻底了解column generation（列生成）算法的原理附java代码
• 列生成和分支定价
• 线性规划技巧: 列生成(Column Generation)
• 列生成算法原理(单纯形基础)
• 列生成（Column Generation）算法
• 【Column Generation思考-02】|从对偶的角度理解Cutting Stock Problem【更新版本】
• Gilmore P C, Gomory R E. A linear programming approach to the cutting-stock problem[J]. Operations research, 1961, 9(6): 849-859.
• 运筹说 第21期 | 算法介绍之列生成算法

前言

之前一直想跟大家分享一下列生成（Column generation），也全网搜了许多文档、视频、论文等。大部分教程抽象程度较高，需要具备大量的基础知识才能看明白，于是写一篇尽可能0基础上手的分享，希望能帮到也在从事相关行业的你。

一定有人会问，有没有行生成算法？当然有啦！！！有机会之后给大家分享！！！——@小猪快跑

从一个例子出发：Cutting Stock Problem

问题描述

卖家有3种长度的木材：9cm（5元），14cm（9元），16cm（10元）。
买家需要木材：4cm（30根），5cm（20根），7cm（40根）。
于是卖家需要通过切割木材满足买家的需求，而且卖家希望成本最低从而达到受益最大。

分析

首先我们来简单思考下，一根9cm的木材能有多少种切法满足买家需要的4cm，5cm，7cm

cost（元）	切前木材长度（cm）	4cm木材数量	5cm木材数量	7cm木材数量
5	9	1	0	0
5	9	2	0	0
5	9	0	1	0
5	9	1	1	0
5	9	0	0	1

本着不浪费的中华传统美德，显然第二行的切法比第一行更好一些哟(＾Ｕ＾)ノ~ＹＯ
我们一般习惯称一种切割方法为cutting pattern。
好啦，之后就非常easy的搞定之后的（只保留比较好的pattern）。

All Cutting Pattern：

cost（元）	切前木材长度（cm）	4cm木材数量	5cm木材数量	7cm木材数量	符号
5	9	2	0	0	$x_1$
5	9	1	1	0	$x_2$
5	9	0	0	1	$x_3$
9	14	3	0	0	$x_4$
9	14	2	1	0	$x_5$
9	14	1	2	0	$x_6$
9	14	1	0	1	$x_7$
9	14	0	1	1	$x_8$
9	14	0	0	2	$x_9$
10	16	4	0	0	$x_{10}$
10	16	2	0	1	$x_{11}$
10	16	1	1	1	$x_{12}$
10	16	0	3	0	$x_{13}$

于是我们只需要从All Cutting Pattern里面确定每个Pattern需要几次就行了。（举个例子：比如${\mathbf{x}}_1=30$代表我们使用第1种切法切了30根9cm长度的木材）。

我们目标是用最少的钱完成买家的需求，也就是每个Pattern的$\mathbf{cost}\times \text{数量}$的总和。

$5x_1+5x_2+5x_3+9x_4+9x_5+9x_6+9x_7+9x_8+9x_9+10x_{10}+10x_{11}+10x_{12}+10x_{13}$

另外呢我们需要满足买家数量上的要求。
4cm要超过30根：

$2x_1+x_2+0x_3+3x_4+2x_5+x_6+x_7+0x_8+0x_9+4x_{10}+2x_{11}+x_{12}+0x_{13}\ge 30$

5cm要超过20根：

$0x_1+x_2+0x_3+0x_4+x_5+2x_6+0x_7+x_8+0x_9+0x_{10}+0x_{11}+x_{12}+3x_{13}\ge 20$

7cm要超过70根：

$0x_1+0x_2+x_3+0x_4+0x_5+0x_6+x_7+x_8+2x_9+0x_{10}+x_{11}+x_{12}+0x_{13}\ge 40$

建模

Master Problem（MP）

根据上面的分析，我们已经有了决策变量$x_i$，我们有了目标函数最少的钱，我们也有了相应的约束满足买家数量上的要求，于是乎我们可以来建模啦！！！

定义决策变量${\mathbf{x}}_{\mathbf{i}}$为使用第$i$种Pattern的根数。
$\min 5x_1+5x_2+5x_3+9x_4+9x_5+9x_6+9x_7+9x_8+9x_9+10x_{10}+10x_{11}+10x_{12}+10x_{13}$
$s.t.$
$2x_1+x_2+0x_3+3x_4+2x_5+x_6+x_7+0x_8+0x_9+4x_{10}+2x_{11}+x_{12}+0x_{13}\ge 30$
$0x_1+x_2+0x_3+0x_4+x_5+2x_6+0x_7+x_8+0x_9+0x_{10}+0x_{11}+x_{12}+3x_{13}\ge 20$
$0x_1+0x_2+x_3+0x_4+0x_5+0x_6+x_7+x_8+2x_9+0x_{10}+x_{11}+x_{12}+0x_{13}\ge 40$
$x_i\in \mathbb{N}$
每一列： 一种可行的切割方案需要被执行的次数，或者多少根棒材使用这种切割方案
每一行：一种棒材的需求必须要被满足

Restricted Master Problem（RMP）

有时候MP的规模会非常大，于是我们考虑如果先用几种Pattern得到一个可行的方案（可能不是最优解），然后我们有了这个可行解后再去依次看如果增加其他Pattern会不会有更好的解，而直接求解MP可能在给定时间内无法得到可行解。

于是顺着这个思路我们构建RMP（要注意我们选择的Pattern至少能得到一个可行解）：

定义决策变量${\mathbf{x}}_{\mathbf{i}}$为使用第$i$种Pattern的根数。
$\min 5x_1+5x_2+5x_3$
$s.t.$
$2x_1+x_2+0x_3\ge 30$
$0x_1+x_2+0x_3\ge 20$
$0x_1+0x_2+x_3\ge 40$
$x_i\in \mathbb{N}$

Restricted Linear Master Problem（RLMP）

常规的列生成方法都是解决LP问题的，这时候我们做一些适当的松弛（也就是$x_i\in \mathbb{N}$变成$x_i\ge 0$）。

于是顺着这个思路我们构建RLMP（要注意RLMP的可行解并不一定是RMP的可行解，但我们只需要把解上取整即可）：

定义决策变量${\mathbf{x}}_{\mathbf{i}}$为使用第$i$种Pattern的根数。
$\min 5x_1+5x_2+5x_3$
$s.t.$
$2x_1+x_2+0x_3\ge 30$
$0x_1+x_2+0x_3\ge 20$
$0x_1+0x_2+x_3\ge 40$
$x_i\ge 0$

然后我们可以用LP求解器直接求解上述模型（这里我用Clp来举例，当然你可以使用你喜欢的求解器）

（为了排版整洁方便阅读，我会把代码放在文末）

直接运行得到

理论上来说我们应该对结果上取整，但我们将在得到松弛后的MP解之后再上取整会更方便高效，所以这里我们不进行上取整（不过好像刚好解是整数Σ(⊙▽⊙"a——@小猪快跑）：
$Objectivevalue=325$
$x_1=5$
$x_2=20$
$x_3=40$
我们试图来解释一下solution的含义：用5个Pattern 1（也就是说选了5根9cm的木材，每根都切成2个4cm），类似的又用了20个Pattern 2和40个Pattern 3，这样的选择能使我们总花费最少——325元。

Dual of Restricted Linear Master Problem

现在我们从另一个角度来考虑RMP。
【卖家】直接购买成品（切割后的木材），转手卖给客户，类似中间商。
【卖家】目标：最多花多少钱购买成品挣得要比自己生产多。
【卖家】约束：需要给成品定一个心理价位。这些心理价位必须满足：对于任意切割方案，购买成品的钱 <= 自己生产的钱。（比如说，卖家花5元买到1根9cm的木材可以切割成2个4cm成品，也就是说卖家如果直接购买2个4cm的成品价格超过5元那不如自己生产了。写成公式：$2y_1\le 5$，其中${\mathbf{y}}_1$为4cm木材的心理价格）
————————————————————————————————————————————————————————
定义决策变量${\mathbf{y}}_{\mathbf{i}}$分别为4cm（30根），5cm（20根），7cm（40根）木材的心理价格。
$\max 30y_1+20y_2+40y_3$
$s.t.$
$2y_1+0y_2+0y_3\le 5$
$y_1+y_2+0y_3\le 5$
$0y_1+0y_2+y_3\le 5$
$y_i\ge 0$

直接求解：

$Objectivevalue=325$
$y_1=2.5$
$y_2=2.5$
$y_3=5.0$
我们试图来解释一下solution的含义：4cm、5cm、7cm成品木材的心理价格分别是2.5元、2.5元、5元，最多花325元购买成品挣得要比自己生产多。也就是说只要价格比2.5元、2.5元、5元低，卖家就赚了。

Subproblem

• 单纯形法检验数的含义：该产品（变量）的市场价格与该产品的隐含成本之差。市场价格高于隐含成本，即检验数大于零时，则可将该产品投入生产。
• 原问题检验数对应其对偶问题的一个基解

如果你知道检验数的概念，会很快理解subproblem，不过在这里我会假设你不知道，试图用最简单的方法让你理解——@小猪快跑

在dual问题中我们计算出成品的心理价格。接下来我们想知道是不是存在更省钱的Pattern。
首先我们来思考一下如何判定一个Pattern是不是更省钱。
例子1：9cm的木材完全浪费掉也是一种Pattern，但一根9cm的木材要5元，也就是说卖家损失了5元。
例子2：Pattern 9：14cm木材（9元）切割成2根7cm成品（使用Pattern 1、2、3的7cm心理价格是5元）。也就是卖家只用Pattern 1、2、3这三种切割方法的时候，只要用低于5元的价格购买7cm再卖掉就比自己切割赚钱了，而如果加上了Pattern 9，卖家如果用5元的价格就不如自己直接用Pattern 9来的赚钱，而是必须低于4.5元（9元/2）的价格购买7cm再卖掉才比自己切割赚钱。
有了上面两个例子我们已经逐渐清晰如何去判定一个Pattern是不是更省钱。
也就是说我们的目标是去寻找一个新Pattern：新Pattern切割后的成品在老Pattern下的心理价格 > 新Pattern的价格
定义决策变量${\mathbf{a}}_{\mathbf{j}}$是新Pattern切割后4cm，5cm，7cm木材的数量。
新Pattern切割后的成品在老Pattern下的心理价格【$2.5a_1+2.5a_2+5a_3$】>新Pattern的价格【9元（14cm）】
即：
$2.5a_1+2.5a_2+5a_3>9$
$\Rightarrow 9-(2.5a_1+2.5a_2+5a_3)<0$（不等号左边我们一般称之为检验数，这就是我们常说的检验数判定解是否是最优解）
这里我们用一个约束转换成目标函数的小技巧：令$z=9-(2.5a_1+2.5a_2+5a_3)$，如果$\min z$比0还大，那说明不存在$z<0$的解，也就是不存在更省钱的Pattern。
新Pattern切割后总长度不能超过原始木材长度：$4a_1+5a_2+7a_3\le 14$
综合起来就是：
$\min 9-(2.5a_1+2.5a_2+5a_3)$
$4a_1+5a_2+7a_3\le 14$

于是我们分别建立9cm、14cm、16cm的Subproblem ：
————————————————————————————————————————————————————————
Subproblem 1（9cm，5元）
定义决策变量${\mathbf{a}}_{\mathbf{j}}$是Pattern切割后4cm，5cm，7cm木材的数量。
$\min z=5-(2.5a_1+2.5a_2+5a_3)$
$s.t.$
$4a_1+5a_2+7a_3\le 9$
$a_i\in \mathbb{N}$
求解得：$a=(2,0,0{)}^T,z=0$
————————————————————————————————————————————————————————
Subproblem 2（14cm，9元）
定义决策变量${\mathbf{a}}_{\mathbf{j}}$是Pattern切割后4cm，5cm，7cm木材的数量。
$\min z=9-(2.5a_1+2.5a_2+5a_3)$
$s.t.$
$4a_1+5a_2+7a_3\le 14$
$a_i\in \mathbb{N}$
求解得：$\mathbf{a}=(0,0,2{)}^{\mathbf{T}},\mathbf{z}=-1$
————————————————————————————————————————————————————————
Subproblem 3（16cm，10元）
定义决策变量${\mathbf{a}}_{\mathbf{j}}$是Pattern切割后4cm，5cm，7cm木材的数量。
$\min z=10-(2.5a_1+2.5a_2+5a_3)$
$s.t.$
$4a_1+5a_2+7a_3\le 16$
$a_i\in \mathbb{N}$
求解得：$a=(4,0,0{)}^T,z=0$

迭代

我们发现Subproblem 2的$z<0$，于是我们把$a=(0,0,2{)}^T$这个Pattern加入到RLMP中，因为这就是我们找到的更省钱的Pattern。
定义决策变量${\mathbf{x}}_{\mathbf{i}}$为使用第$i$种Pattern的根数。
$\min 5x_1+5x_2+5x_3+9x_9$
$s.t.$
$2x_1+x_2+0x_3+0x_9\ge 30$
$0x_1+x_2+0x_3+0x_9\ge 20$
$0x_1+0x_2+x_3+2x_9\ge 40$
$x_i\ge 0$

求解后得到：

$Objectivevalue=305$
$x_1=5$
$x_2=20$
$x_3=0$
$x_9=20$
$y_1=2.5$
$y_2=2.5$
$y_9=4.5$
由于$x_3=0$，所以我们可以从模型中舍去$x_3$。
于是我们分别建立9cm、14cm、16cm的Subproblem ：
————————————————————————————————————————————————————————
Subproblem 1（9cm，5元）
定义决策变量${\mathbf{a}}_{\mathbf{j}}$是Pattern切割后4cm，5cm，7cm木材的数量。
$\min z=5-(2.5a_1+2.5a_2+4.5a_9)$
$s.t.$
$4a_1+5a_2+7a_9\le 9$
$a_i\in \mathbb{N}$
求解得：$a=(2,0,0{)}^T,z=0$
————————————————————————————————————————————————————————
Subproblem 2（14cm，9元）
定义决策变量${\mathbf{a}}_{\mathbf{j}}$是Pattern切割后4cm，5cm，7cm木材的数量。
$\min z=9-(2.5a_1+2.5a_2+4.5a_9)$
$s.t.$
$4a_1+5a_2+7a_9\le 14$
$a_i\in \mathbb{N}$
求解得：$a=(0,0,2{)}^T,z=0$
————————————————————————————————————————————————————————
Subproblem 3（16cm，10元）
定义决策变量${\mathbf{a}}_{\mathbf{j}}$是Pattern切割后4cm，5cm，7cm木材的数量。
$\min z=10-(2.5a_1+2.5a_2+4.5a_9)$
$s.t.$
$4a_1+5a_2+7a_9\le 16$
$a_i\in \mathbb{N}$
求解得：$a=(4,0,0{)}^T,z=0$

所有的Subproblem目标$z\ge 0$，那就说明不存在更好的Pattern，于是最终的解（记得上取整）：$x_1=5,x_2=20,x_9=20$

以上整个算法过程我们就称之为——列生成

列生成：Cutting Stock Problem

本节需要一定的运筹学基础，但如果你已经看完上文的话，我相信理解起来也会非常简单了。——@小猪快跑

提出Gomory割的大佬Gomory在IBM的时候，与另一个大佬Gilmore 共同提出了著名的Column Generation：

参考文献：
@article{gilmore1961linear,
title={A linear programming approach to the cutting-stock problem},
author={Gilmore, Paul C and Gomory, Ralph E},
journal={Operations research},
volume={9},
number={6},
pages={849–859},
year={1961},
publisher={INFORMS}
}

问题描述

卖家有$n$种长度的木材：$L_i$ m（$c_i$元）
买家需要$m$种长度的木材：$l_i$ m（$d_i$根）
于是卖家需要通过切割木材满足买家的需求，而且卖家希望成本最低从而达到受益最大。

建模

Master Problem（MP）

$$\begin{array}{rll}\min & \sum _{i=1}^n{c}_i{x}_i& \\ s.t.& \sum _{i=1}^n{a}_{ij}{x}_i⩾b_i,& j=1,2,...,m\\ & x_i\in \mathbb{N},& \forall i\end{array}$$

Restricted Master Problem（RMP）

$$\begin{array}{rll}\min & \sum _{i=1}^k{c}_i{x}_i& \\ s.t.& \sum _{i=1}^k{a}_{ij}{x}_i⩾b_i,& j=1,2,...,m\\ & x_i\in \mathbb{N},& \forall i\end{array}$$

也就是MP问题变量数从$n$减少到$k$个，需要注意我们强制了$x_i(i>k)$的变量为非基变量。

Dual of Restricted Master Problem

为了在Subproblem中计算检验数：${\sigma }_j=c_j-c_B{B}^{-1}{a}_j$，我们需要计算出$c_B{B}^{-1}$，一般来说他有两重含义：

1. 通过求解RMP问题得到的影子价格（shadow price）。
2. 通过求解RMP对偶问题得到的对偶变量(dual variable)。

我们一般使用第2种对偶问题求解。因为RMP一般变量极多，而单纯形法对于变量多的问题求解很困难。原问题的变量对应对偶问题的约束，目标系数对应约束边界，约束矩阵倒转过来。所以在对偶问题上就没有这个困扰。

Subproblem

$$\begin{array}{rl}\min & c_j-\sum {\omega }_i{a}_i\\ s.t.& \sum _{i=1}^m{l}_i{a}_i⩽L_j\\ & a_i⩾0,\forall i\end{array}$$

迭代

检验数小于0的话，就把Pattern加入列，进行下一次迭代，直到所有检验数大于等于0。最后上取整结果即可。

流程图

总结

1. 列生成算法主要用于求解变量多，但大部分变量取值为0的线性规划问题。
2. 总体思路是先选出小部分变量构建RMP（也就是其余变量取值都是0）快速得到一个可行解，之后再通过Subproblem计算检验数（reduced cost）寻找有没有更好的变量加入模型再次求解，直到找不到更好的变量。
3. 为什么增加变量可以让目标值更好呢？原问题增加变量=>对偶问题增加约束=>对偶问题最优值不变或者变小（对偶问题是max问题，约束越多可行域越小，自然目标函数优度下降）=>原问题最优值不变或者变小（对偶问题和原问题最优解是一样的）
   事实上，我们需要找的变量是能在对偶问题中让最优解不满足新增的约束。

代码

from gurobipy import GRB, Model, quicksum, Column

EPS = 1e-6


def solve_cutting_stock_with_column_generation(widths_to_demand, material_width_to_price):
    """
    使用列生成方法解决切割库存问题。
    :param widths_to_demand: 需求宽度及其数量的字典
    :param material_width_to_price: 原材料宽度及其价格的字典
    :return: 找到的所有模式
    """
    widths = list(widths_to_demand.keys())
    demand = list(widths_to_demand.values())

    # 初始化模型
    m = Model("Restricted Linear Master Problem")
    m.setAttr(GRB.Attr.ModelSense, GRB.MINIMIZE)
    # 取消模型打印信息
    m.setParam("OutputFlag", 0)

    # 主问题初始化
    init_material_width = min(material_width_to_price)
    patterns_material_width = [init_material_width] * len(widths)
    patterns_price = [material_width_to_price[init_material_width]] * len(widths)
    patterns = []  # 存储所有已知的模式
    for i, w in enumerate(widths):
        tmp = [0] * len(widths)
        tmp[i] = init_material_width // w
        patterns.append(tmp)

    # 添加变量
    var_x = [
        m.addVar(
            lb=0,
            ub=GRB.INFINITY,
            obj=price,
            vtype=GRB.CONTINUOUS,
            name=f"x_{pattern}".replace(" ", "")
        ) for price, pattern in zip(patterns_price, patterns)
    ]

    # 添加约束条件
    cons = [
        m.addConstr(
            quicksum(pattern[i] * x for x, pattern in zip(var_x, patterns)) >= d,
            name=f"cons_{i}".replace(" ", "")
        ) for i, d in enumerate(demand)
    ]

    epoch = 1
    while True:
        print("-" * 20, f"Iteration {epoch}", "-" * 20)
        # 重置模型状态
        m.reset()

        # 求解主问题
        m.optimize()
        print(f"RLMP Obj: {m.ObjVal}")
        print(f"RLMP Var: {m.getAttr('X', m.getVars())}")

        # 如果主问题不可行或已经找到最优解，退出循环
        if m.Status != GRB.OPTIMAL:
            if m.Status == GRB.INFEASIBLE:
                print("主问题无可行解")
            elif m.Status == GRB.UNBOUNDED:
                print("主问题无界")
            else:
                print(f"主问题状态: {m.Status}")
            break

        # 删除不需要的pattern变量
        # drop_idx = [i for i, x in enumerate(var_x) if x.X < EPS]
        # for i in drop_idx[::-1]:
        #     patterns.pop(i)
        #     patterns_material_width.pop(i)
        #     patterns_price.pop(i)
        #     m.remove(var_x.pop(i))

        # 获取对偶变量
        duals = m.getAttr("Pi", m.getConstrs())
        print(f"Dual Values: {duals}")

        # 求解子问题来找到新的模式
        flag = True
        for material_width, price in material_width_to_price.items():
            new_pattern = solve_subproblem(widths, material_width, price, duals)
            # 如果找到了新模式，添加到列表中
            if new_pattern is not None:
                patterns.append(new_pattern)
                patterns_material_width.append(material_width)
                patterns_price.append(material_width_to_price[material_width])
                var_x.append(m.addVar(
                    lb=0,
                    ub=GRB.INFINITY,
                    obj=price,
                    vtype=GRB.CONTINUOUS,
                    column=Column(new_pattern, cons),
                    name=f"x_{new_pattern}".replace(" ", "")
                ))
                print(f"New Pattern: {new_pattern}")
                flag = False
        if flag:
            break

        epoch += 1

    # 将RMP转化为MIP
    print("-" * 20, f"Solve MIP", "-" * 20)
    m.setAttr("VType", m.getVars(), GRB.INTEGER)
    m.optimize()
    print(f"MIP Obj: {m.ObjVal}")
    print(f"MIP Var: {m.getAttr('X', m.getVars())}")

    return patterns


def solve_subproblem(widths, material_width, price, duals):
    # 创建子问题模型
    m = Model("subproblem")
    m.setParam("OutputFlag", 0)
    var_x = [m.addVar(vtype=GRB.INTEGER, name=f"x_{i}") for i in range(len(widths))]

    # 目标函数：最大化剩余材料的价值
    m.setObjective(price - quicksum(d * x for d, x in zip(duals, var_x)), GRB.MINIMIZE)

    # 约束条件：总宽度不超过材料宽度
    m.addConstr(quicksum(w * x for w, x in zip(widths, var_x)) <= material_width, "width_constraint")

    m.optimize()

    if m.Status == GRB.OPTIMAL and m.objVal < 0:
        pattern = [int(round(v.X)) for v in var_x]
        return pattern
    else:
        if m.Status != GRB.OPTIMAL:
            print(f"子问题状态: {m.Status}")
        return None


if __name__ == '__main__':
    # 示例数据 默认原材料宽度 > 需求宽度
    material_width_to_price = {9: 5, 14: 9, 16: 10}  # 原材料宽度 - 价格
    widths_to_demand = {4: 30, 5: 20, 7: 39}  # 需求宽度 - 需求数量

    # 调用函数
    patterns = solve_cutting_stock_with_column_generation(widths_to_demand, material_width_to_price)
    print("Patterns:", patterns)

-------------------- Iteration 1 --------------------
RLMP Obj: 370.0
RLMP Var: [15.0, 20.0, 39.0]
Dual Values: [2.5, 5.0, 5.0]
New Pattern: [1, 1, 0]
New Pattern: [1, 2, 0]
New Pattern: [0, 3, 0]
-------------------- Iteration 2 --------------------
RLMP Obj: 320.0
RLMP Var: [5.0, 0.0, 39.0, 20.0, 0.0, 0.0]
Dual Values: [2.5, 2.5, 5.0]
New Pattern: [0, 0, 2]
-------------------- Iteration 3 --------------------
RLMP Obj: 300.5
RLMP Var: [5.0, 0.0, 0.0, 20.0, 0.0, 0.0, 19.5]
Dual Values: [2.5, 2.5, 4.5]
-------------------- Solve MIP --------------------
MIP Obj: 301.0
MIP Var: [5.0, -0.0, 1.0, 20.0, -0.0, -0.0, 19.0]
Patterns: [[2, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0], [1, 2, 0], [0, 3, 0], [0, 0, 2]]