列生成算法及python实践

列生成算法及python实践

列生成算法是在混合整数规划中一个经典的且强大的精确算法。

引例-下料问题

一家木材厂有一种长度为80英寸的木材原料。现在有3个顾客，他们分别需要46个10英寸，22个11英寸，43个19英寸的木材产品。
为了满足客户不同的需求并且使得切割的原料最少，问该公司最少需要的木材多少根，以及如何切割？

为了方便求解该问题，我们引入如下记号：

记号	含义
$m$	客户的集合
$W$	木材的最大宽度
$w_i，i\in m$	客户$i$要求的棒材宽度
$n_i，i\in m$	客户$i$对棒材宽度为$w_i$的需求数量

求解方法

Kantorovich 建模方式

记号	含义
$\mathcal{K}$	待切割的棒材的集合
y k ∈ { 0 , 1 } y_k\in\{0,1\} yk​∈{0,1}	棒材$k$是否被切割，1表示切割，0表示不切割
$x_i^k$	棒材$k$被规格为$w_i$切割的个数

于是可以将数学模型表述为如下形式：
min ∑ k ∈ K y k s.t. ∑ k ∈ K x i k ≥ n i , i = 1 , … , m , (demand) ∑ i = 1 n w i x i k ≤ W y k , k ∈ K , (width limitation) x i k ∈ Z + , y k ∈ { 0 , 1 } . \begin{aligned} &\text{min} \sum_{k\in\mathcal{K}}y_k \\ &\text{s.t.} \sum_{k\in\mathcal{K}}x_i^k\geq n_i,\quad i=1,\ldots,m,\quad\text{(demand)} \\ &\sum_{i=1}^nw_ix_i^k\leq Wy_k,\quad k\in\mathcal{K},\quad\text{(width limitation)} \\ &x_i^k\in\mathbb{Z}_+, y_k\in\{0,1\}. \end{aligned} ​mink∈K∑​yk​s.t.k∈K∑​xik​≥ni​,i=1,…,m,(demand)i=1∑n​wi​xik​≤Wyk​,k∈K,(width limitation)xik​∈Z+​,yk​∈{0,1}.​

缺点

• 线性松弛界质量太差，当问题规模大时需要构造大量的二进制变量和整数变量。

列生成算法

列生成的建模方式与Kantorovich建模方式不同，需要引入pattern的概念，即对于一根木材原料切割方案，比如一根80英寸的原料可以切割成8个10英寸的产品。下表给出了一些pattern：

	10英寸	11英寸	19英寸	浪费的木材长度	切割的次数
pattern1	8	0	0	0	$x_1$
pattern2	6	1	0	9	$x_2$
pattern3	6	0	1	1	$x_3$
pattern4	5	1	1	0	$x_4$
客户的需求	46	22	43	0

pattern的数量会随着需求产品种类的增加而指数级增加，针对上表中的pattern，我们给出4种切割方案的模型：
$$\begin{array}{rl}\text{min}& \sum _{i=1}^{4}1\ast x_i\\ \text{s.t.}& 8x_1+6x_2+6x_3+5x_4\ge 46\\ & 0x_1+1x_2+0x_3+1x_4\ge 22\\ & 0x_1+0x_2+1x_3+0x_4\ge 43\\ & x_1,x_2,x_3,x_4\in {\mathbb{Z}}_{+}\end{array}$$
其实上述的模型并没有把pattern给枚举完全，也即模型的列不完全，所以被称为限制性主问题（restrictive master problem）。单纯形法的经验告诉我们求解像上述的3行的优化问题的基是一个3×3的矩阵。故我们需要寻找到"最有用"的列作为基，而单纯形法中判断一列是不是有用，就是判断该列的检验数是否是负的，针对最小化问题。如果某列的检验数是负的，入基则会减少目标函数的值。所有的非基变量的检验数都大于等于0的话，此时的解为最优解。我们把寻找新列的问题叫子问题或者定价问题。针对$j$列的检验数为：
$$reducecost=c_j-{\mathbf{c}}_{\mathbf{B}}^{\mathbf{T}}{\mathbf{B}}^{-1}{\mathbf{N}}_{\mathbf{j}}$$
其中${\mathbf{c}}_{\mathbf{B}}^{\mathbf{T}}{\mathbf{B}}^{-1}$是原问题的对偶变量的值，$N_j$是$j$列的基向量，在本问题中代表某种切割方案。设原问题的对偶变量为$\mathbf{y}$,切割方案为${\mathbf{N}}_{\mathbf{j}}=(a_1,a_2,a_3{)}^{\top }$. 那么子问题可以写作：
$$\begin{array}{rl}\text{min}& 1-{\mathbf{y}}^{\mathbf{T}}(a_1,a_2,a_3{)}^{\top }\\ \text{s.t.}& 10a_1+11a_2+19a_3\le 80\\ & a_1,a_2,a_3\in {\mathbb{Z}}_{+}\end{array}$$
列生成算法的步骤为：

• 初始化：为主问题生成一个矩阵A，可以使用启发式的方法；
• 重复：
• 1. 求解限制性线性松弛的主问题，令${\mathbf{y}}^{\top }={\mathbf{c}}_{\mathbf{B}}^{\mathbf{T}}{\mathbf{B}}^{-1}$为该问题的对偶变量值。
• 2. 求解背包问题类型的子问题，其最优值为$k$。
• 3. 将子问题的最优解变成新列的系数加入限制性线性松弛的主问题。
• 直到$k\ge 0$

列生成的缺点：

• 线性松弛的最优解不一定原问题的最优解。

实验结果

问题规模	求解方式	最优目标函数值	时间s
n=3	kantorovich	20	0.2093
n=3	column generation	20	0.01
n=10	kantorovich	61	3.97
n=10	column generation	63	0.01

总结

1. kantorovich的建模方式求解速度较慢，但求解结果更准确。
2. column generation的建模方式求解速度快，但求解结果不一定准确，获得是全局最优解的上界，介意这个gap的话可以采用分支定价算法。

本文python代码地址：https://github.com/gaopan812/OR_practice/tree/main/%E5%88%97%E7%94%9F%E6%88%90%E7%AE%97%E6%B3%95

参考文献

[1]孙小玲,李端.整数规划[M].科学出版社,2010.
[2]运筹优化常用模型、算法及案例实战——Python+Java实现