列生成算法简介

1. 什么是列生成

列生成算法是一种用于解决大规模线性规划问题的高效算法，它基于单纯形法的思想，通过求解子问题来找到可以进基的非基变量。在列生成算法中，每个变量都代表一列，因此称为列生成算法。该算法的优点在于其高效的计算性能和较好的收敛性，适用于处理大规模、复杂的线性规划问题。

在列生成算法的迭代过程中，因为会不断有变量入基，所以会导致限制主问题的列不断增加，所以叫做列生成算法。

2. 列生成的应用范围

列生成被广泛应用于调度问题、切割问题、车辆路径问题、选址问题等。 该算法的优点在于其高效的计算性能和较好的收敛性，适用于处理大规模、复杂的线性规划问题。对于变量数目很多的线性优化问题，单纯形法速度很慢，可以用到列生成方法来加快求解速度。

3. 列生成的原理

基本思路如下：

1、先把原问题限制到一个规模更小的限制主问题，在限制主问题的基础上用单纯形法求解，但此时的解并不是主问题的最优解。

2、通过一个子问题去检查那些未被考虑的变量中是否有使得reduced cost小于0的？如果有，就把这个变量的相关系数列加入到限制主问题的系数矩阵中，回到第一步。

经过反复迭代，知道子问题的reduced cost rate大于等于0，那么主问题就求到了最优解。

4. 基本概念

受限主问题

$$
min(y_1+y_2+\cdots+y_k)
$$

$s.t.$
$$
R_1: a_{11}y_1+\cdots+a_{1k}y_k \geq b_1
$$
$$
R_2: a_{21}y_1+\cdots+a_{2k}y_k \geq b_2
$$

$$
\cdots
$$

$$
R_m:a_{m1}y_1+\cdots+a_{mk}y_k \geq b_m
$$

就是从主问题中选取k个变量构成的松弛问题。

子问题

通过求解RMP问题或者RMP对偶问题后，得到了想要的$c_BB^{-1}$以后，子问题就是通过$\sigma_j = c_j - c_BB^{-1}a_j$

这条公式，在$y_{k+1}、y_m$中寻找检验数为负且最小的变量，将变量对应的那一列添加到RMP中。

列生成算法实现

案例求解

问题描述

cutting stock problem是采用列生成算法求解的经典案例。该问题如下：有一些纸筒，每个纸筒长度为16m。顾客需要25个3m，20个6m，15个7m的纸筒。要求在满足顾客需求的情况下，裁剪的纸筒数最小。

建立模型

可以采用启发式算法求解到一个初始解如下：

5个纸筒采用切割方案1：3，3，3，3，3；

10个纸筒采用切割方案2：6，6；

8个纸筒采用切割方案3：7，7；

总计23个纸筒。可以看出，采用23个纸筒是一定可以满足要求的，这是问题的一个上界（upper bound）。

该问题可行的切割方案很多，我们可以采用P表示所有可行裁剪方案的集合，里面的方案总数为n（这里的n并不需要知道其确切取值，一般来说n是一个很大的数）。$a_{ij}$表示第$j$种方案里类别$i$的个数，$y_j$表示第$j$种方案的选择个数。建立数学模型如下：
$$
min \quad y_1+y_2+\cdots+y_n
$$
$s.t.$
$$
R1:a_{11}y_1 + \cdots + a_{1n}y_n \geq 25
$$
$$
R2:a_{21}y_1 + \cdots + a_{2n}y_n \geq 20
$$

$$
R3:a_{31}y_1 + \cdots + a_{3n}y_n \geq 15
$$

约束中的每一列对应的是一种切割方案，前三种方案已知，就是我们前面采用启发式算法求解的三种方案。其它的$n-3$种切割方案未知。

问题求解

第一轮循环

首先从上述模型中选出一些列，构成问题的限制主问题，这里我们选取前3列。则限制主问题如下
$$
min \quad y_1 + y_2 + y_3
$$
$$
5y_1 + 0 y_2 + 0y_3 \geq 25
$$

$$
0y_1 + 2 y_2 + 0y_3 \geq 20
$$

$$
0y_1 + 0 y_2 + 2y_3 \geq 15
$$

使用python+gurobi求解上述线性规划问题

import gurobipy as gp
# 实例化模型
m = gp.Model('column generation')
y = m.addVars(3,vtype = gp.GRB.CONTINUOUS,name='y',lb=0)
m.addConstr(5*y[0] >= 25)
m.addConstr(2*y[1] >= 20)
m.addConstr(2*y[2] >= 15)
m.setObjective(y[0]+y[1]+y[2],gp.GRB.MINIMIZE)
m.setParam('OutputFlag',0)
m.optimize()
m.getAttr(gp.GRB.Attr.Pi)

可以求解出对偶变量的取值为$c_BB^{-1}=[0.2,0.5,0.5]$。现在要找一列加入RMP，但是并不知道其取值，我们将其记作$\alpha_4 = [a_{14},a_{24},a_{34}]^T$。得到非基变量的检验数为$\sigma_4 = c_4 - c_BB^{-1}\alpha_4 = 1-0.2\alpha_{14} - 0.5\alpha_{24} - 0.5\alpha_{34}$

从而构造子问题如下：
$$
min(1-0.2a_{14} - 0.5a_{24} - 0.5a_{34})
$$
$s.t.$
$$
3a_{14}+6a_{24}+7a_{34} \leq 16
$$
$$
a_{ij} \in Z
$$

同样采用python+gurobi求解上述问题

import gurobipy as gp
m1 = gp.Model('subproblem')
a = m1.addVars(3,vtype=gp.GRB.INTEGER,lb=0,name='a')
m1.setObjective(1-0.2*a[0]-0.5*a[1]-0.5*a[2])
m1.addConstr(3*a[0]+6*a[1]+7*a[2] <= 16)  # 列生成规则
m1.setParam('OutputFlag',0)
m1.optimize()
print(a[0].x,a[1].x,a[2].x)

求解出$\alpha_4 = [1,2,0]^T$，检验数为$\sigma_4 = c_4 - c_BB^{-1}\alpha_4 = 1-0.2\times1 - 0.5\times2 - 0.5\times0 = -0.2 < 0$，因为检验数小于0，所以需要将$y_4$入基，添加到主问题中，开始第二轮迭代。

第二轮循环

加入$y_4$后，限制主问题变为
$$
min \quad y_1 + y_2 + y_3 + y_4
$$
$$
5y_1 + 0 y_2 + 0y_3 + 1y_4\geq 25
$$

$$
0y_1 + 2 y_2 + 0y_3 + 2y_4 \geq 20
$$

$$
0y_1 + 0 y_2 + 2y_3 + 0y_4 \geq 15
$$

限制主问题与之前相比，多了一列，所以被称为列生成算法。求解该限制主问题，其对偶变量取值为$c_BB^{-1}=[0.2,0.4,0.5]$，下一个需要入基的变量记为$\alpha_5 = [a_{15},a_{25},a_{35}]^T$，构造子问题如下
$$
min(1-0.2a_{15} - 0.4a_{25} - 0.5a_{35})
$$
$s.t.$
$$
3a_{15}+6a_{25}+7a_{35} \leq 16
$$
$$
a_{ij} \in Z
$$

采用同样方法求解该子问题，得到结果$\alpha_5 = [3,0,1]^T$，检验数为$\sigma_5 = c_5 - c_BB^{-1}\alpha_5 = 1-0.2\times3 - 0.4\times0 - 0.5\times1 = -0.1 < 0$，因为检验数小于0，所以需要将$y_5$入基，添加到主问题中，开始下一轮迭代。

第三轮循环

加入$y_5$后，限制主问题变为
$$
min \quad y_1 + y_2 + y_3 + y_4 + y_5
$$
$$
5y_1 + 0 y_2 + 0y_3 + 1y_4 + 3y_5\geq 25
$$

$$
0y_1 + 2 y_2 + 0y_3 + 2y_4 + 0y_5 \geq 20
$$

$$
0y_1 + 0 y_2 + 2y_3 + 0y_4 + 1y_5 \geq 15
$$

还是求解限制主问题，其对偶变量取值为$c_BB^{-1}=[0.1667,0.4167,0.5]$，下一个需要入基的变量记为$\alpha_6 = [a_{16},a_{26},a_{36}]^T$，构造子问题如下
$$
min(1-0.1667a_{16} - 0.4167a_{26} - 0.5a_{36})
$$
$s.t.$
$$
3a_{16}+6a_{26}+7a_{36} \leq 16
$$
$$
a_{ij} \in Z
$$

求得结果$\alpha_5 = [1,1,1]^T$，检验数$\sigma_6= c_6 - c_BB^{-1}\alpha_6 = 1-0.1667\times1- 0.4167\times1 - 0.5\times1 = -0.08333 < 0$，所以将$y_6$代入

第四轮循环

加入$y_6$后，限制主问题变为
$$
min \quad y_1 + y_2 + y_3 + y_4 + y_5 + y_6
$$
$$
5y_1 + 0 y_2 + 0y_3 + 1y_4 + 3y_5 + y_6\geq 25
$$

$$
0y_1 + 2 y_2 + 0y_3 + 2y_4 + 0y_5 + y_6 \geq 20
$$

$$
0y_1 + 0 y_2 + 2y_3 + 0y_4 + 1y_5 + y_6 \geq 15
$$

还是求解限制主问题，其对偶变量取值为$c_BB^{-1}=[0.2,0.4,0.4]$，下一个需要入基的变量记为$\alpha_7 = [a_{17},a_{27},a_{37}]^T$，构造子问题如下
$$
min(1-0.2a_{17} - 0.4a_{27} - 0.4a_{37})
$$
$s.t.$
$$
3a_{17}+6a_{27}+7a_{37} \leq 16
$$
$$
a_{ij} \in Z
$$

求得结果$\alpha_5 = [5,0,0]^T$，检验数$\sigma_6= c_6 - c_BB^{-1}\alpha_6 = 1-0.2\times5- 0.4\times0 - 0.4\times0 = 0$，结束迭代，此时列生成算法结束。

此时，我们将$y_7$代入模型，此时求解结果就是最优解。

求解出最优解为$y=[1,0,0,3,1,14,0]$，此时目标函数值为19。

所以我们得到的最终切割方案为

方案	该方案的数量
3，3，3，3，3，3	1
3，6，6	3
3，3，3，7	1
3，6，7	14

完整版列生成代码

import gurobipy as grb
from gurobipy import GRB
import numpy as np


def master_problem(column, vtype):
    m = grb.Model()
    x = m.addMVar(shape=column.shape[1], lb=0, vtype=vtype)
    m.addConstr(lhs=column @ x >= demand_number_array)
    m.setObjective(x.sum(), GRB.MINIMIZE)
    m.optimize()

    if vtype == GRB.CONTINUOUS:
        return np.array(m.getAttr('Pi', m.getConstrs()))
    else:
        return m.objVal, np.array(m.getAttr('X'))


def restricted_lp_master_problem(column):
    return master_problem(column, GRB.CONTINUOUS)


def restricted_ip_master_problem(column):
    return master_problem(column, GRB.INTEGER)


def knapsack_subproblem(kk):
    m = grb.Model()
    x = m.addMVar(shape=kk.shape[0], lb=0, vtype=GRB.INTEGER)
    m.addConstr(lhs=demand_width_array @ x <= roll_width)
    m.setObjective(1 - kk @ x, GRB.MINIMIZE)
    m.optimize()

    flag_new_column = m.objVal < 0
    if flag_new_column:
        new_column = m.getAttr('X')
    else:
        new_column = None
    return flag_new_column, new_column


roll_width = np.array(16)
demand_width_array = np.array([3, 6, 7])
demand_number_array = np.array([25, 20, 15])
initial_cut_pattern = np.diag(np.floor(roll_width / demand_width_array))


flag_new_cut_pattern = True
new_cut_pattern = None
cut_pattern = initial_cut_pattern
while flag_new_cut_pattern:
    if new_cut_pattern:
        cut_pattern = np.column_stack((cut_pattern, new_cut_pattern))
    kk = restricted_lp_master_problem(cut_pattern)
    flag_new_cut_pattern, new_cut_pattern = knapsack_subproblem(kk)

minimal_stock, optimal_number = restricted_ip_master_problem(cut_pattern)

print('************************************************')
print('parameter:')
print(f'roll_width: {roll_width}')
print(f'demand_width_array: {demand_width_array}')
print(f'demand_number_array: {demand_number_array}')
print('result:')
print(f'minimal_stock: {minimal_stock}')
print(f'cut_pattern: {cut_pattern}')
print(f'optimal_number: {optimal_number}')

运行结果

************************************************
roll_width:16
demand_width_array:[3,6,7]
demand_number_array:[25,20,15]
minimal_stock:19.0
cut_pattern:[[5,0,0,1,3,1]
[0,2,0,2,0,1]
[0,0,2,0,1,1]]
optimal_number:[1,0,0,3,1,14]