线性规划技巧之列生成

在制造业中，常常要切割原材料，比如把纸卷、布料、钢管等等，切割成不同长度，从而满足不同产品，或者不同用户的需求。

由于需求长度不同，切割方式有多种。不同的切割方式，消耗的原料总数不同。问题是怎么切割，既能满足需求，又最节省原料。这样的问题，称为 下料问题 (The cutting stock problem)。本文以下料问题为例介绍列生成算法，它可以用来解这类线性规划问题。

例子

已知原材料长度为 1000。有四种需求类型，其需求长度和需求量如下表所示。

现在要回答两个问题：一是最少需要多少根原材料；二是每根原材料怎么切割。

为了方便描述，先定义切割方式。用向量表示，如下所示：

$$ a_j = \begin{bmatrix} a_{1,j}\\ a_{2,j}\\ a_{3,j}\\ a_{4,j} \end{bmatrix} $$

其中 $j$ 是切割方式的编号，$a_{ij}$ 代表第 $i$ 种长度的切割数量，$i=1,2,3,4$。

例如，$(2, 0, 0, 0)$ 表示按第一种需求（长度为450）切两个；$(0, 1, 0, 4)$ 表示按第二种需求（长度360）切一个，再按第四种需求（长度140）切四个。

需要注意的是，还要保证切割是“可行的”。比如说，$(3, 0, 0, 0)$ 不可行，因为总长度只有1000，切不成三段长度450。

具体来说，令 $L$ 代表原材料长度，向量 $s$ 代表需求长度，即 $$ s = \begin{bmatrix} s_1\\ s_2\\ s_3\\ s_4 \end{bmatrix} $$ 满足如下条件： $$ a_j^T s \leq L $$

可以这样理解，一根原材料对应一个切割向量。用 $x_j$ 代表按 $a_j$ 切割的原材料数量，那么 $\sum_j x_j$ 就是所需要的原材料总数。

回到问题，它可以这样描述：已知所有切割方式 $a_j$，要计算它的个数 $x_j$，使得需求被满足，且所需原材料总数最小。

说明一下，所有的切割方式（满足可行性），可以通过枚举得到。

数学模型

有了前面的分析，不难构建出数学模型。

指标

• $i$ — 需求
• $j$ — 切割方式

参数

• $d_i$ — 需求量
• $a_{i,j}$ — 在第 $j$ 种切割方式中，需求长度 $i$ 的切割数量

决策变量

• $x_j \in \mathbb{N}$ — 第 $j$ 种切割方式需要的原材料数量

目标

• $\min \sum_j x_j$ — 最小化原材料数量

约束

• 满足需求量：$\sum_{j} a_{i,j} x_j \geq d_i$

数学规划模型 $$ \begin{aligned} \min ~ & \sum_j x_j \\ \text{ s.t. } & \sum_{j} a_{i,j} x_j \geq d_i, ~\forall i \\ & x_j \in \mathbb{N}, ~\forall j \end{aligned} $$

注意，上式是整数规划，本身就不好解。即使是它的线性规划，也不一定好解。原因是它的变量非常多，比约束多得多。

下面简单分析一下。用 $m$ 代表需求种类数，容易发现，它等于约束数量。在这个例子中，$m$ 等于 4。

设 $n$ 是变量个数，它等于切割方式的数量。需求长度有 $O(2^m)$ 种组合，如果都是可行的，那么切割方式的数量 $n \in O(2^m)$。

这意味着，问题规模随 $m$ 呈指数增长，那么求解时间也是这样。

求解思路

考虑到整数规划比较难，尤其是大规模问题，最优解很难求。下面的思路是求近似解。

整个过程分为四步：

第一步，考虑整数规划的松弛形式，即把 $x_j\in \mathbb{N}$ 松弛为 $x_j\geq 0$，即 $$ \begin{aligned} \min ~ & \sum_j x_j \\ \text{ s.t. } & \sum_{j} a_{i,j} x_j \geq d_i, ~\forall i \\ & x_j \geq 0, ~\forall j \end{aligned} $$

第二步，求解上述线性规划，得到松弛解 $x^$。这里有个问题，$x^$ 可能不是整向量，也就是说不可行（相对原问题）。

说明：当问题规模较大时，直接求解困难。下面介绍列生成方法，可以用来解这样的问题。

第三步，把松弛解取整，即 $\bar{x}= \lfloor x^* \rfloor$。但是还有个问题，取整后可能导致需求约束不满足，此时依然不可行。

第四步，满足剩下的需求，得到可行解。比如说，按需求的长度从大到小排序，然后依次满足。

第三和第四步比较简单，这里不多做解释。下面介绍第二步，用列生成求解上述线性规划。

列生成

回顾单纯形算法，它从一个基本解出发，朝着目标函数值下降的方向迭代。

现在的问题是，变量规模太大。或者说，系数矩阵的列太多，导致需要的内存太大，或者计算时间过长。

但它有个特点，很多可行的切割方式，其实不在最优解中。换句话说，大部分的列，可能不需要考虑。

列生成的思想是，只考虑必要的列。初始化的时候，从一组基开始，迭代一次增加一列。新增的列就是新的入基，通过求解子问题得到。

下面定义主问题和子问题，其中主问题用来算可行解，子问题用来生成列。

主问题

原问题可以写成矩阵形式： $$ \begin{aligned}\min~ & c^T x \ & Ax = b \ & x \geq 0 \end{aligned} $$ 其中 $A$ 是系数矩阵，每一列表示一种可行的切割方式，$c$ 和 $b$ 的定义如下： $$ c= \begin{bmatrix}1\\ 1\\ \vdots\\ 1\end{bmatrix}, ~ b = \begin{bmatrix}d_1\\ d_2\\ \vdots \\ d_m \end{bmatrix} $$ 主问题 的形式不变，区别是初始化的时候，系数矩阵是任意的基矩阵 $B$，即 $A:=B$。换句话说，主问题一开始只考虑 $m$ 列。

在这个例子中，可以令 $$ B =\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} $$ 容易验证，$B$ 的每一列都是可行的切割方式。

子问题

子问题的任务是，生成主问题系数矩阵的列。要满足两个条件：第一是可行性，即生成的列是可行的切割方式；第二是有效性，即生成的列作为入基，要使目标函数值下降。

具体来说，令 $y\in\mathbb{N}^m$ 代表子问题的解，也就是要生成的列。

先看可行性，要满足如下条件： $$ s^T y \leq L $$ 其中 $s\in\mathbb{R}^m $ 是每种需求的长度，$L\in\mathbb{R}$ 是原材料的长度。

再看有效性，目标函数下降的方向，就是梯度小于零的方向。回顾单纯形算法，目标函数的梯度也称作 Reduced Cost，即 $$ \mu^T := c^T_N - \lambda^T N $$ 其中 $\lambda^T = c^T_B B^{-1}$ 是对偶变量值，也称 Shadow Price；$B$ 和 $N$ 对应基本矩阵和非基矩阵；$c_B$ 和 $c_N$ 作为系数，对应基变量和非基变量。

为了满足有效性，一个自然的想法是，让 $y$ 对应的梯度分量 $\mu_y$ 尽可能小，即 $$ \min~ \mu_y := c_y - \lambda^T y $$ 当 $\mu_y < 0$ 时，说明目标函数还可以下降，因此满足有效性，否则原问题达到最优。

在下料问题中 $c = (1, 1, …, 1)^T$，我们有 $$ \mu_y= 1 -\lambda^T y $$ 这样一来， $$ \min~ 1-\lambda^T y \Leftrightarrow \max~ \lambda^Ty $$ 结合上述两个条件，我们得到 子问题： $$ \begin{aligned} \max~ & \lambda^T y \\ \text{s.t. } & s^T y \leq L \\ & y \in \mathbb{N}^m \end{aligned} $$ 说明一下，主问题求解之后，$\lambda$ 已经计算出来，无需重复计算。

迭代过程

根据上面的例子，我们介绍算法的流程。

第一步，初始化主问题，其中系数矩阵 $A$ 考虑一组基： $$ A :=\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} $$

第二步，求解主问题： $$ \begin{aligned}\min~ & c^Tx \\ & Ax = b \\ & x \geq 0 \end{aligned} $$ 其中 $$ c= \begin{bmatrix}1\\ 1\\ \vdots \\ 1\end{bmatrix}, ~ b = \begin{bmatrix}97 \\ 610 \\ 395 \\ 211 \end{bmatrix} $$ 注意，$c$ 和 $x$ 的维度与系数矩阵 $A$ 列的数量相同。当 $A$ 增加列时，$c$ 和 $x$ 也要相应增加。

第三步，求解子问题： $$ \begin{aligned} \max~ & \lambda^T y \\ \text{s.t. } & s^T y \leq L \\ & y \in \mathbb{N}^m \end{aligned} $$ 其中 $s^T=(450,360,310,140)$，代表四种需求长度；$\lambda$ 来自主问题的对偶变量值。

第四步，判断是否最优。

如果 $\mu_y = 1-\lambda^T y \geq 0$，那么达到最优，停止迭代并输出最优解。

否则，把 $y$ 添加到主问题的列，即 $$ A:= \begin{bmatrix} A & y \end{bmatrix} $$ 然后回到第二步（继续迭代）。

总结

本文以下料问题为例，介绍了列生成方法。它是一种基于单纯形算法的间接求解方法，适合这种指数多个变量的问题。相比直接求解，列生成可以降低问题规模，让问题变得可解。求解的过程，就是不断求解主问题和子问题，直到满足停止条件。最后回到下料问题，解线性规划只是第二步，还要把它变成整数解，以及满足需求量约束，从而得到可行解。