【Column Generation思考-02】|从对偶的角度理解Cutting Stock Problem【更新版本】

本文致谢：感谢【运小筹读者2群】的小伙伴们与我就该话题展开的热烈而有用的讨论。以及幸新杰老师提供的有用的信息。

Cutting Stock Problem是运筹优化中一个非常重要的问题。该问题与一个非常强大的算法框架【Column Generation】直接相关。本节我们从对偶的角度，来理解Cutting Stock Problem。

另外，很多人懂的Column Generation的过程，但是也许有些读者并不了解Column Generation迭代过程中，问题的Lower bound如何计算。因为一些读者在Column Generation的过程中，仅仅关注子问题的reduced cost是否为负，不需要去关心每一步Column Generation的迭代中，全局的Upper bound和Lower bound。因为在CG结束迭代的时候，CG reformulation得到的全局Lower bound恰好等于当前RMP的线性松弛RLMP的LP relaxation解。也就是此时$$\text{global LB}=z_{\text{RLMP}}$$
所以大家不关心也没关系。不过，知道Lower bound如何计算，会让你对整个算法理解更上一层楼，这篇推文，我们就与大家共同推导一下Column Generation求解Cutting Stock Problem的过程中，Lower bound如何计算。

Cutting Stock Problem简介

参考文献：

@article{Berberler2011,
title = {A software for the one-dimensional cutting stock problem},
journal = {Journal of King Saud University - Science},
volume = {23},
number = {1},
pages = {69-76},
year = {2011},
issn = {1018-3647},
doi = {https://doi.org/10.1016/j.jksus.2010.06.009},
url = {https://www.sciencedirect.com/science/article/pii/S1018364710000741},
author = {Murat Erşen Berberler and Urfat Nuriyev and Ahmet Yıldırım},
}

下料问题

一家企业拥有一批原料棒材，他要将这些棒材进行切割加工成为规定尺寸的产品，然后卖给顾客。这家企业需要设计最优的棒材切割方案，以最小的原材料消耗，生产出规定尺寸的产品，以满足顾客的需求。

如上图（例子来自文献Berberler et al., 2011），一家公司有长度为$L$的原料棒材。其客户需要的产品尺寸和需求量分别为

产品规格	需求量
$\frac{2L}{3}$	1
$\frac{L}{2}$	1
$\frac{L}{3}$	4
$\frac{L}{4}$	2

下面我们给出该问题的一般模型。首先引入下面的参数和决策变量：

参数

• $L$: 原料棒材长度；
• $\mathcal{K}$: 可用棒材的集合，以$k$索引；
• $\mathcal{I}$: 需求产品的集合，$i\in \mathcal{I}$表示第$i$种需求的index；
• $d_i$: 第$i$种需求的需求量；
• $h_i$: 第$i$种产品的长度。

决策变量

• $y_k$: 第$k$根棒材是否被切割；
• $x_{ik}$: 第$k$根棒材使用的切割方案中，包含的第$i$种产品的数量；

基于此，Cutting Stock Problem可以被建模为下面的整数规划模型：

min ⁡ ∑ k ∈ K y k s . t . ∑ k ∈ K x i k ⩾ d i , ∀ i ∈ I , ∑ i ∈ I h i x i k ⩽ L y k , ∀ k ∈ K , x i k , y k ∈ { 0 , 1 } , ∀ i ∈ I , ∀ k ∈ K . \begin{aligned} \min \quad & \sum_{k \in \mathcal{K}} y_k && \\ s.t. \quad & \sum_{k \in \mathcal{K}} x_{ik} \geqslant d_i, && \forall i \in \mathcal{I}, \\ & \sum_{i \in \mathcal{I}} h_ix_{ik} \leqslant Ly_k, && \forall k \in \mathcal{K}, \\ & x_{ik}, y_k \in \{0, 1\} , && \forall i \in \mathcal{I}, \forall k \in \mathcal{K}. \end{aligned} mins.t.​k∈K∑​yk​k∈K∑​xik​⩾di​,i∈I∑​hi​xik​⩽Lyk​,xik​,yk​∈{0,1},​​∀i∈I,∀k∈K,∀i∈I,∀k∈K.​
其中

• 约束1： 表示所有切割得到的产品，必须满足每一种需求；
• 约束2： 每一根棒材切割产生的产品的长度，必须小于单根棒材的总长度。这个是保证切割方案可行性。如果棒材$k$不被使用，则其不会得到任何产品。

Cutting Stock Problem的Column generation reformulation

Cutting Stock Problem的上述模型求解效率并不会特别高。但是可以通过Column generation达到求解效率的显著提升。提出Gomory割的大佬Gomory在IBM的时候，与另一个大佬Gilmore 共同提出了著名的Column Generation（参加下面文献）。我们目前不做特别详细的描述，直接给出Column Generation求解Cutting Stock Problem的模型。

参考文献：
@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}
}

主问题 （Master Problem, MP）

假设所有的可行切割方案均已知，我们令其为$P$，引入下面的参数和决策变量

参数

• $a_{ip}$: 第$p$种切割方案中，第$i$种需求产品的数量。

决策变量

• $x_p$: 执行第$p$种切割方案的棒材的数量。

$$\begin{array}{rlrl}\min & \sum _{p\in P}{x}_p& & \\ s.t.& \sum _{p\in P}{a}_{ip}{x}_p⩾d_i,& & \forall i\in \mathcal{I},\\ & x_p⩾0\text{ and integer},& & \forall p\in P.\end{array}$$

当需求产品的个数较少时，也就是$|\mathcal{I}|$较小时，$P$可以比较容易被穷举出来。但是$|\mathcal{I}|$较大，$L$较大时，穷举出$P$中的所有元素，就不再简单。因此我们可以使用Column Generation的方法，找出那些比较promising的切割方案，加入到主问题中，摒弃那些unpromising的切割方案。promising指的是，将这个切割方案加入到主问题中，会使得主问题的目标函数下降。在Column Generation中，就是这个切割方案的reduced cost小于0。

因此，在Column Generation迭代中，我们只是考虑了一部分已经得到的切割方案。由于一个切割方案对应一个决策变量，也就是对应约束矩阵中的一列，因此一个切割方案又叫做一列。我们的主问题，由于只考虑了$P$的一个子集，我们用$P^{\prime }$表示CG过程中主问题包含的列。基于$P^{\prime }$的主问题模型，又被称之为Restricted Master Problem（RMP，限制性主问题）。为了得到对偶变量，我们将其松弛，松弛后的RMP我们称其为RLMP。RLMP的数学模型如下：
$$\begin{array}{rlrl}\min & \sum _{p\in P^{\prime }}{x}_p& & \\ s.t.& \sum _{p\in P^{\prime }}{a}_{ip}{x}_p⩾d_i,& & \forall i\in \mathcal{I},\\ & x_p⩾0,& & \forall p\in P^{\prime }.\end{array}$$

我们用$\pi$表示RLMP的约束的对偶变量。

Column Generation的子问题就是找到reduced cost为负的列。这一过程也叫pricing。Column Generation求解Cutting Stock Problem的子问题就是找到一个切割方案，并且其reduced cost为负。我们这里可以尝试找到具有最负检验数的方案。因此模型可以写为

$$\begin{array}{rlrl}\min & 1-\sum _{i\in \mathcal{I}}{a}_i{\pi }_i& & \\ s.t.& \sum _{i\in \mathcal{I}}{h}_i{a}_i⩽L& & \\ & a_i⩾0\text{ and integer},& & \forall i\in \mathcal{I}.\end{array}$$

Cutting Stock Problem的小例子+Python调用Gurobi求解

我们考虑下面的小例子：

棒材	Value
棒材	20
产品长度	9和6
需求量	20个6英寸的产品和 15个9英寸的产品

6: 20
9: 15 

我们穷举出所有可行的切割方案，建立等价的MP。

$$\begin{array}{r}\begin{array}{ccccccccccccc}\min & x_1& +& x_2& +& x_3& +& x_4& +& x_5& +& x_6& \\ 6英寸& 3x_1& +& 0x_2& +& x_3& +& 2x_4& +& 1x_5& +& 0x_6& ⩾20\\ 9英寸& 0x_1& +& 2x_2& +& x_3& +& 0x_4& +& 0x_5& +& 1x_6& ⩾15\\ & x_1,& & x_2,& & x_3& +& x_3& +& x_3& +& x_3& ⩾0\text{ and integer}\end{array}\end{array}$$

将其松弛成线性规划：

$$\begin{array}{r}\begin{array}{ccccccccccccc}\min & x_1& +& x_2& +& x_3& +& x_4& +& x_5& +& x_6& \\ 6英寸& 3x_1& +& 0x_2& +& x_3& +& 2x_4& +& 1x_5& +& 0x_6& ⩾20\\ 9英寸& 0x_1& +& 2x_2& +& x_3& +& 0x_4& +& 0x_5& +& 1x_6& ⩾15\\ & x_1,& & x_2,& & x_3& +& x_3& +& x_3& +& x_3& ⩾0\end{array}\end{array}$$

并且调用Gurobi求解LMP（MP的线性松弛）：

from gurobipy import * 
model = Model('CG')
x = {}
for i in range(1, 7):
    x[i] = model.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name = 'x_' + str(i))
    
model.setObjective(x[1] + x[2] + x[3] + x[4] + x[5] + x[6], GRB.MINIMIZE)

model.addConstr(3*x[1] + 0*x[2] + x[3] + 2*x[4] + 1*x[5] + 0*x[6] >= 20)
model.addConstr(0*x[1] + 2*x[2] + x[3] + 0*x[4] + 0*x[5] + 1*x[6] >= 15) 

model.optimize()

for key in x.keys():
    print('x[{}] = {}'.format(key, x[key].x)) 

运行结果

Solved in 0 iterations and 0.01 seconds
Optimal objective  1.416666667e+01
x[1] = 6.666666666666667
x[2] = 7.5
x[3] = 0.0
x[4] = 0.0
x[5] = 0.0
x[6] = 0.0

可见线性松弛的解为14.17，是个小数解。

我们再将其设置为整数，直接求解MP本身，结果为

x[i] = model.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name = 'x_' + str(i))

Out:[]
Optimal solution found (tolerance 1.00e-04)
Best objective 1.500000000000e+01, best bound 1.500000000000e+01, gap 0.0000%
x[1] = 7.0
x[2] = 8.0
x[3] = -0.0
x[4] = -0.0
x[5] = -0.0
x[6] = -0.0

可以看到，下界14.1，上界15，所以15就是最优解了，因为最优解一定是一个整数。

因此，这里就可以得到一个结论，也许很多小伙伴并不了解。

列生成的formulation: 如果列出了所有的方案(列), 求解最终的LMP一定会得到整数解吗？

• 答案： 不一定。即使完整地写出了MP，最后求解LMP，得到的解可能是小数，也可能是整数。
  但是，如果是小数，这个RLMP提供的下界也是个与最优解差别在1以内的小数。比如，LMP的解是$14.16$，但是最优解一定是整数，那就最少是$15$。结果我们把决策变量设置成整数以后，得到的解正好是$15$。因此，也就相当于Gap = 0了。

另外，我们也可以看看最终的LMP和其对偶问题的解是否相同：

• 原问题：
  $$\begin{array}{r}\begin{array}{ccccccccccccc}\min & x_1& +& x_2& +& x_3& +& x_4& +& x_5& +& x_6& \\ 6英寸& 3x_1& +& 0x_2& +& x_3& +& 2x_4& +& 1x_5& +& 0x_6& ⩾20\\ 9英寸& 0x_1& +& 2x_2& +& x_3& +& 0x_4& +& 0x_5& +& 1x_6& ⩾15\\ & x_1,& & x_2,& & x_3& +& x_3& +& x_3& +& x_3& ⩾0\end{array}\end{array}$$
  每一列： 一种可行的切割方案
  每一行：一种棒材的需求必须要被满足

• 对偶问题
  $$\begin{array}{r}\begin{array}{ccccc}\max & 20y_1& +& 15y_2& \\ & 3y_1& +& 0y_2& ⩽1\\ & 0y_1& +& 2y_2& ⩽1\\ & 1y_1& +& 1y_2& ⩽1\\ & 2y_1& +& 0y_2& ⩽1\\ & 1y_1& +& 0y_2& ⩽1\\ & 0y_1& +& 1y_2& ⩽1\\ & y_1,& & y_2& ⩾0\end{array}\end{array}$$

我们用Python调用Gurobi求解LMP的对偶，即求解Dual LMP，代码如下：

from gurobipy import * 
model = Model('Dual LMP')
y = {}
for i in range(1, 3):
    y[i] = model.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name = 'y_' + str(i))
    
model.setObjective(20 * y[1] + 15 * y[2], GRB.MAXIMIZE)

model.addConstr(3*y[1] + 0*y[2] <= 1)
model.addConstr(0*y[1] + 2*y[2] <= 1)
model.addConstr(1*y[1] + 1*y[2] <= 1)
model.addConstr(2*y[1] + 0*y[2] <= 1)
model.addConstr(1*y[1] + 0*y[2] <= 1)
model.addConstr(0*y[1] + 1*y[2] <= 1)

model.optimize()

for key in y.keys():
    print('y[{}] = {}'.format(key, y[key].x)) 

设置为小数，求解结果如下

Solved in 0 iterations and 0.01 seconds
Optimal objective  1.416666667e+01
y[1] = 0.3333333333333333
y[2] = 0.5

设置为整数，求解结果如下

Optimal solution found (tolerance 1.00e-04)
Best objective -0.000000000000e+00, best bound -0.000000000000e+00, gap 0.0000%
y[1] = 0.0
y[2] = 0.0

课件，Dual LMP在设置成整数的时候，最优解居然为0。

因此这里我们可以稍微做一点讨论，我们考虑下面四个模型。

Primal IP
$$\begin{array}{rl}\min & cx\\ s.t.& Ax⩾b\\ & x⩾0\text{ and integer}\end{array}$$

将IP 松弛 得到Primal LP
Primal LP
$$\begin{array}{rl}\min & cx\\ s.t.& Ax⩾b\\ & x⩾0\end{array}$$

将Primal LP对偶，得到Dual LP。
Dual LP
$$\begin{array}{rl}\max & by\\ s.t.& yA⩽c\\ & y⩾0\end{array}$$

将Dual LP设置为整数规划，得到Dual IP。
Dual IP
$$\begin{array}{rl}\max & by\\ s.t.& yA⩽c\\ & y⩾0\text{ and integer}\end{array}$$

我们可以得到以下关系：

$$z_{\text{Dual IP}}⩽z_{\text{Dual LP}}⩽z_{\text{Primal LP}}⩽z_{\text{Primal IP}}$$

在强对偶成立的时候
$$z_{\text{Dual LP}}=z_{\text{Primal LP}}$$

但是，我们想知道$z_{\text{Dual IP}}$和$z_{\text{Primal IP}}$的差距究竟有多大？
根据上述例子，我们可知
$$\begin{array}{rl}& z_{\text{Primal IP}}=15\\ & z_{\text{Dual IP}}=0\end{array}$$
差距也太大了。所以$z_{\text{Dual IP}}$和$z_{\text{Primal IP}}$的关系，还是真的不好评判。

好了，这个小分享就到此为止。

我们来进入今天的重要内容：如何直观解释Cutting Stock Problem的对偶问题？

其实笔者一直认为，每一个原问题的对偶问题，都应该有比较合理的直观解释。它会向我们提供个其他视角的建模理念或者思路。

直观解释CG求解Cutting Stock Problem对偶问题的含义

CG求解Cutting Stock Problem的原问题

模型的含义含义：

• 原问题：
  $$\begin{array}{r}\begin{array}{ccccccccccccc}\min & x_1& +& x_2& +& x_3& +& x_4& +& x_5& +& x_6& \\ 6英寸& 3x_1& +& 0x_2& +& x_3& +& 2x_4& +& 1x_5& +& 0x_6& ⩾20\\ 9英寸& 0x_1& +& 2x_2& +& x_3& +& 0x_4& +& 0x_5& +& 1x_6& ⩾15\\ & x_1,& & x_2,& & x_3& +& x_3& +& x_3& +& x_3& ⩾0\end{array}\end{array}$$
  每一列： 一种可行的切割方案需要被执行的次数，或者多少根棒材使用这种切割方案
  每一行：一种棒材的需求必须要被满足

原问题直观的含义就是：

比如这个生产厂家就是【A公司】

• 原问题：A公司的客户下了一批订单，要6英寸的产品20个，9英寸的棒材15个。A公司需要向外面购买棒材，来生产这些需求。其中，每一根棒材的成本为1美元。当然，问题里其实是假设他的棒材已经买好了的。不过我们其实可以看做是他还没买呢，他先做个计划，然后决策出来最少的原材料成本之后，再去购买卷钢。
• 原问题的目标函数，就是最小化购买卷钢的成本，它等于$\min {\sum }_{r\in R}{c}_r\cdot x_r$，只不过这里每一根卷钢的购买成本都是一样的，都是1美元。并且满足购买的卷钢可以生产出规定量需求的产品。当然了，$\min {\sum }_{r\in R}\cdot x_r$就代表买的用于执行各个切割方案的钢材的总成本。

• 其对偶问题的思考角度，就可以从外购的视角来重新思考这个问题。首先我们给出对偶问题的形式（这里直接以案例为例给出对偶问题）。

• 对偶问题的模型如下：
  $$\begin{array}{r}\begin{array}{ccccc}\max & 20y_1& +& 15y_2& \\ & 3y_1& +& 0y_2& ⩽1\\ & 0y_1& +& 2y_2& ⩽1\\ & 1y_1& +& 1y_2& ⩽1\\ & 2y_1& +& 0y_2& ⩽1\\ & 1y_1& +& 0y_2& ⩽1\\ & 0y_1& +& 1y_2& ⩽1\\ & y_1,& & y_2& ⩾0\end{array}\end{array}$$

从外购的视角，上述对偶问题直观的含义就可以解释为：

对偶问题仍然可以看做是【A公司做决策】

• 【A公司】觉得，公司可以将这些生产任务可以外包，公司不自己生产。如果外包的话，公司需要付给承担外包生产任务的企业相应的加工费，这个加工费就是生产一个6英寸产品的单价和生产一个9英寸产品的价格(其实就是心理价格)。或者可以理解为：公司就直接买现成的产品，即：向其他公司购买现成的产品，然后转手交付给客户。比如，你向我买，我再去向别人买，我就是个中间商，哈哈。
• 如果这个【A公司】想做个中间商的角色，他想买现成的产品交差，他需要决策：我最多可以允许多少的预算，能够使得我买现成的交差比我自己生产要划算。所以他要$\max$自己最多允许花多少钱，因为这样他好掂量掂量哪种方式更划算。其中，$y_1$和$y_2$就是他给自己定的6英寸和9英寸产品的心理价位。这些心理价位必须满足：
  每一种切割方案下，我购买现成的产品花的钱，比直接买原材料棒材自己生产要节省。
  如果买现成的产品花的钱比我自己直接买卷钢自己生产花的多，那我不如自己买卷钢自己生产了。
• 所以就有
  $3\cdot y_1+0\cdot y_2⩽1$
  这个约束。其余约束的解释类似。

综上所述，对偶问题的目标函数和决策变量的直观含义如下：

• 目标函数： $\max 20\cdot y_1+15\cdot y_2$，表示，最大化自己直接购买现成产品交差(或者外包给其他方生产)的总成本，其中，每单位产品的价格，是按照心理价位来计算。也就是外包/外购的时候，一单位的6英寸的产品的价格$y_1$以及一单位的9英寸的产品的价格$y_2$。之所以$\max$，是因为想找到一个自己的最大开销的边界，找到这个边界之后，我就可以决定，现成的产品价格在多少以内，我去直接买现成的产品更划算。现成的产品价格超过什么范围之后，我就自己生产更划算。
• 约束：$3\cdot y_1+0\cdot y_2⩽1$，表示：我外购现成的产品花费的钱，要比我自己买卷钢生产这些产品的成本要小。否则，我不如自己购买卷钢，自己生产。其中：$3\cdot y_1+0\cdot y_2$就是外购现成的产品花费的钱；$1$就是自己买卷钢自己生产的成本。如果我要决定外购，我必须要保证，外购对我来讲更省钱。
• 比如上述代码，找到了对偶问题$\max$的最优解$$\begin{gathered} y_1=0.333; \\ {y}_2=0.5 \end{gathered}$$，目标函数值$$Z=14.1667$$。这说明，我们外购的最大预算为$Z=14.1667$，如果外购现成产品花的钱超过$Z=14.1667$，则说明我们就亏了，我们不如买卷钢自己生产。并且对应的，临界的6英寸的产品和9英寸的产品的外购心理单位价格为$y_1=0.333;y_2=0.5$，如果超过其他生产商的现成的产品的价格超过这个数字，我们就不外购，我们就自己买卷钢自己生产更划算。如果说现成的产品价格低于这个数字，那说明我们可以直接购买现成的产品进行交付，因为直接外购更划算。

其他问题的对偶视角解释

桌椅生产问题

原问题：

• 目标函数： 最大化总收入，得到每种产品的生产量
• 约束1： 木材的使用量<=总的可用量： 生产桌子用的木材+生产椅子用的木材<=总的木材可用量
• 约束2： 钢材的使用量<=总的可用量： 生产桌子用的钢材+生产椅子用的钢材<=总的钢材可用量

原问题：

• 目标函数： 最小的出租收入：决策每种原材料的出租价格
• 约束1： 生产一张桌子的原材料卖的总钱数>=自己生产赚的钱
• 约束2： 生产一把椅子的原材料卖的总钱数>=自己生产赚的钱

总结：生产问题

• 原问题变量：生产多少
• 原问题约束：可用资源上限
• 对偶问题变量：出租价格
• 对偶题约束：出租收入大于等于自己生产的收入

最短路问题

另外一个例子就是最短路问题(具体模型推导本公众号之前的推文)。

总结：最短路问题

• 原问题变量：路径$x_{ij}$
• 原问题约束：路径连贯性约束
• 对偶问题变量：到达每个点的最短距离${\pi }_i$
• 对偶题约束：一条弧上的终点的最短距离-起点的最短距离<= d[i,j]

参考文献

• [1]. Berberler M E, Nuriyev U, Yıldırım A. A software for the one-dimensional cutting stock problem[J]. Journal of King Saud University-Science, 2011, 23(1): 69-76.
• [2]. Gilmore P C, Gomory R E. A linear programming approach to the cutting-stock problem[J]. Operations research, 1961, 9(6): 849-859.
• [3]. 刘兴禄, 熊望祺, 臧永森, 段宏达, 曾文佳, 陈伟坚. 运筹优化常用模型、算法及案例实战：Python+Java实现. 清华大学出版社（即将出版）.