【大规模整数规划】| Benders Decomposition: 一份让你满意的【入门-编程实战-深入理解】的学习笔记

作者：刘兴禄，清华大学，清华-伯克利深圳学院

初次完稿日期：2020.10.08

欢迎关注我们的微信公众号 运小筹

大规模混合整数规划求解算法的综述

大规模混合整数规划(mixed integer programming, MIP)的求解在运筹优化中具有至关重要的地位。这里我们强调一下，不做特殊说明的话，我们考虑的都是混合整数线性规划(mixed integer linear programming, MILP)。求解MIP的常用精确算法主要包括：

• branch and bound
• cutting plane
• branch and cut
• column generation(不保证得到最优解，需要结合其他技巧获得最优解)
• branch and price
• Lagrangian relaxation(不保证得到最优解，需要结合其他技巧获得最优解)
• Dantzig-Wolfe decomposition
• Benders decomposition
• …

其他衍生的拓展方法，本文不做介绍。

其中，branch and bound、branch and cut是非常通用的求解框架，适用于所有的MIP 。

Lagrangian relaxation并不能保证得到最优解，但是运气好的时候，也可以得到最优解。不过Lagrangian relaxation可以得到比线性松弛（LP relaxation）更紧的界限(Bound)，因此可以用来加速branch and bound或者branch and cut。

而column generation和branch and price这两个框架的通用性就弱很多。这两个框架都是基于模型重构的。只有可以重构为相应问题的模型，才可以使用这两种方法。此外，column generation不能保证获得全局最优解（具体原因可以查看相关文献，也可以阅读笔者即将出版的教材，见参考文献1）。而将branch and bound的框架嵌入column generation中，构成branch and price框架，则可以保证获得全局最优解。

以上精确算法中：

• branch and bound是一种分而治之的思想，本质上是一种隐枚举（branching）。但是由于加入了prune(剪枝)和bounding（定界）的步骤，使得该框架在实际求解中要比真正的纯枚举要高效得多。
• branch and cut是目前最流行的求解MIP的求解器的通用算法框架。由于其在branch and bound的基础上，加入了cutting plane的步骤，用于割去当前节点的最优小数解，从而逼近该节点的凸包，从而显著地加速了求解过程。
• column generation和branch and price是基于模型重构而来的算法。通常是将原来的MIP分解成为一个主问题（Master Problem）和若干个子问题（Subproblem），然后迭代求解。当然，子问题又叫定价子问题（Pricing Problem）。由于分解之后，各个部分的求解相对容易，以及主问题一般是更为紧凑的模型，因此会提供更好的界限，从而加速收敛。这些框架都是融合了模型重构和分解的思想。
• Lagrangian relaxation是一种松弛的思想。通过结合对偶理论，得到比单纯的线性松弛更好的界限。
• Dantzig-Wolfe decomposition和Benders decomposition是根据模型的特殊结构，将模型分解为更容易求解的小问题，通过小问题之间的迭代求解和交互，最终达到精确求解模型的目的的精确算法。但是二者的分解思路并不相同。Dantzig-Wolfe decomposition是基于一个表示定理得来的分解方法，该方法需要MIP的约束矩阵符合块角状的特征，通用性有限，使用之前需要考察模型是否符合该结构。而Benders decomposition实际上是较为通用的分解框架，CPLEX中也包含了Benders decomposition的算法组件，用户可以自己定制分解策略。Benders decomposition的主要思想是，将较复杂的模型分解成为2部分，分别求解2部分都较为容易，通过两部分之间的交互迭代，最终算法得到收敛，得到最优解。

每一种精确算法都不存在绝对的优劣关系，它们各有特点，都蕴含着精妙的思想。在实际情况中，我们需要灵活应用，巧妙结合，从而达到显著加速求解的目的。

本文就来着重介绍分解算法中的重要成员: Benders decomposition。这里强调一下，本文涉及的Benders decomposition，是最基本的Benders decomposition，更高级别的Benders decomposition的拓展算法，请读者自行阅读相关文献。

分解算法

上面讲到，分解算法，是基于模型的特点，将大问题拆解称为若干小问题。单独求解这些小问题的难度，远低于直接求解原来的大模型。这些小问题之间需要进行一定的交互，从而保证整分解后的模型能够收敛到原来模型的最优解，从而保证分解方法的全局最优性。

这里需要明确，分解$\ne$拆解。分解的操作，实际上是很有艺术性和理论性的。

分解包括拆解$+$耦合。其中拆解就是将问题拆成若干小问题，耦合就是将各个小问题联系起来，形成交互，从而保证全局最优性。

也许目前的直观解释还不够具体，不过没关系，我们先说清楚整体的思想，然后再以具体的理论加以解释。

Benders Decomposition的算法原理

考虑如下的MIP模型：
min ⁡ x , y c x + f y s.t. A x + B y = b x ⩾ 0 y ∈ Y ⊆ R n

$$\begin{aligned} \underset{x,y}{\min} \quad &cx+fy \\ \text{s.t.} \quad &Ax + By = b \\ & x\geqslant 0 \\ & y \in Y \subseteq \mathbb{R}^n \end{aligned}$$ x,ymin​s.t.​cx+fyAx+By=bx⩾0y∈Y⊆Rn​
其中$c\in {\mathbb{R}}^{1\times m},f\in {\mathbb{R}}^{1\times n}$,是行向量，$x,y$是列向量，且$x,y$是决策变量，其维度分别为$m\times 1$和$n\times 1$。$A,B$是约束矩阵，其维度分别为$A\in {\mathbb{R}}^{r\times m},B\in {\mathbb{R}}^{r\times n}$。$b$为右端常数，其维度为$b\in {\mathbb{R}}^{r\times 1}$。注意这里的符号，${\mathbb{R}}^{r\times m}$其实就是表示是一个$r\times m$的一个实数矩阵。

该MIP模型具有下面的特点：

• $x$是连续型(continuous)决策变量，较为简单；
• $y$是复杂决策变量，可以认为是$0-1$型(binary)决策变量或者整数(integer)型决策变量，相较$x$而言，略复杂，因此在模型中我们用了数学语言$y\in Y\subseteq {\mathbb{R}}^n$来表达。前一个$y\in Y$，这个$Y$可以是一系列的 y ∈ { 0 , 1 } y \in \{0,1\} y∈{ 0,1}的约束，也可以是$y\text{ integer}$的约束。当然，还可以是$y⩾0$的连续约束，都可以。为了方便，我们统一用$y\in Y$来表示。

由于$x$为简单变量，$y$为复杂变量，当规模较大时，直接求解上述MIP比较困难。这个困难的原因是：

• 复杂变量$y$搅合进来了，这个罪魁祸首给问题求解带来了挑战。

因此我们设想：如果我们将这个捣蛋者$y$先搞定，剩下的部分中，只有$x$是决策变量，问题就变成了线性规划(Linear programming, LP)了，这不就简单太多了嘛？

Benders Decomposition就是基于这样的思想。

我们观察到，如果给定$y$的值，假定为$\bar{y}$，那么剩余部分的模型就可以写成
SP: min ⁡ x c x s.t. A x = b − B y ˉ x ⩾ 0

$$\begin{aligned} \text{SP:} \quad \underset{x}{\min} \quad &cx \\ \text{s.t.} \quad &Ax = b - B\bar{y} \\ & x\geqslant 0 \end{aligned}$$ SP:xmin​s.t.​cxAx=b−Byˉ​x⩾0​

该模型是LP，求解难度比之前的MIP降低了好几档次。因为MIP是NP-hard问题，而LP是多项式时间可精确求解的，不是NP-hard。该问题一般被称之为子问题（Subproblem problem, SP）。

如何给出$y$的值$\bar{y}$呢？我们可以将关于$x$的部分全部忽略掉，变成下面的模型
MP 0 : min ⁡ y f y s.t. y ∈ Y ⊆ R n

$$\begin{aligned} \text{MP}_0: \quad \underset{y}{\min} \quad & fy \\ \text{s.t.} \quad & y \in Y \subseteq \mathbb{R}^n \end{aligned}$$ MP0​:ymin​s.t.​fyy∈Y⊆Rn​

求解${\text{MP}}_0$非常容易，求解得到的解，即可作为$\bar{y}$传给$\text{SP}$。
该问题一般被称之为主问题（Masterproblem, MP）。${\text{MP}}_0$加了下标0表示初始的$\text{MP}$(因为后续迭代过程，MP需要加入cut)。

但是，分别求解模型$\text{MP}$和$\text{SP}$并不能等同于求解了原MIP。要想达到等价地求解原MIP的目的，还需要$\text{MP}$和$\text{SP}$之间的交互。

为什么要交互呢？

原因是：

• $\text{MP}$删除了所有关于$x$的约束，相当于抛弃了$x$的影响。但是实际上$x$一定会影响$\text{MP}$。我们通过一种先试错，再补救的过程，先删除所有的$x$的相关部分，构成初始的$\text{MP}$，通过求解$\text{SP}$，获得一些信息，这些信息反映了$x$对$y$的影响，我们通过某种方式，将$x$对$y$的影响加回到$\text{MP}$中，进行补救。通过整个过程的迭代，我们最终就可以达到求解原MIP的目的。

其中，补救的步骤，实际上专业术语叫做recourse，反应在具体形式上，就是以cutting plane的形式，将补救措施加回到$\text{MP}$中。

至于如何补救补救，那就需要用到对偶理论的知识。首先说明，Benders Decomposition分解的补救措施，是以两种cutting plane的形式加回到$\text{MP}$中的，分别为

• Benders optimality cut和
• Benders feasibility cut。

这两种Cut来源于$\text{SP}$的对偶，由于$\text{SP}$是LP，因此我们可以将$\text{SP}$对偶，变成$\text{Dual SP}$，其形式如下：
Dual SP: max ⁡ α α T ( b − B y ˉ ) s.t. A T α ⩽ c α  free

$$\begin{aligned} \text{Dual SP:} \quad \underset{\alpha}{\max} \quad & \alpha^T (b-B\bar{y}) \\ \text{s.t.} \quad &A^T\alpha \leqslant c \\ & \alpha \text{ free} \end{aligned}$$ Dual SP:αmax​s.t.​αT(b−Byˉ​)ATα⩽cα free​
我们可以观察到：

• $\text{Dual SP}$的可行性不依赖于$\bar{y}$的取值。
• 如果强对偶成立(即$\text{SP}$和$\text{Dual SP}$都有可行解)，则$z_{\text{SP}}=z_{\text{Dual SP}}$。
• 根据弱对偶性，有$z_{\text{SP}}⩾z_{\text{Dual SP}}$，即$\text{Dual SP}$为$\text{SP}$提供了一个下界。

这里我们不加证明的给出Benders optimality cut和Benders feasibility cut的具体形式。（为什么不加证明呢？是小编不会嘛？当然不是，是太麻烦了，小编在自己写的教材里面给出了非常详细的论述，大家可以去自行查看。或者可以阅读下面的参考文献。）

• Benders optimality cuts的具体形式为
  ( α p i ) T ( b − B y ) ⩽ q , ∀ i ∈ P

  $$\begin{aligned} (\alpha_p^i)^T(b-By) \leqslant q, \quad \forall i \in \mathcal{P} \end{aligned}$$ (αpi​)T(b−By)⩽q,∀i∈P​
  其中，$\mathcal{P}$表示迭代过程中找到的子问题的对偶问题$\text{Dual SP}$的极点(extreme point，其实就是最优解)的集合，${\alpha }_p^i$就表示$\text{Dual SP}$的一个极点。

• Benders feasibility cuts的具体形式为
  ( α r i ) T ( b − B y ) ⩽ 0 , ∀ i ∈ R

  $$\begin{aligned} (\alpha_r^i)^T(b-By) \leqslant 0, \quad \forall i \in \mathcal{R} \end{aligned}$$ (αri​)T(b−By)⩽0,∀i∈R​
  其中，$\mathcal{R}$表示迭代过程中找到的子问题的对偶问题$\text{Dual SP}$的极射线(extreme ray，其实就是无界射线)的集合，${\alpha }_r^i$就表示$\text{Dual SP}$的一个极射线。

这里，我们省去具体的理论论证部分，直接给出简洁的结论。完整，详细的理论介绍，参见文献

• Zeki Caner Taşkin. Benders Decomposition. American Cancer Society, 2011. ISBN 9780470400531.

此外，笔者即将出版的教材《运筹优化常用模型、算法与案例实战——Python+Java实现》第17章Benders分解算法给出了更为详尽、直观的解释，可以帮助读者更容易理解详细的原因。待出版之后，有兴趣的读者可以去查看。

Benders Decomposition的完整步骤如下：

1. 我们将原问题拆分为一个主问题（Master problem, MP）和一个子问题（Subproblem problem, MP）。$\text{MP}$如下
   $$\begin{array}{rl}\text{MP}:\underset{y}{\min }& fy+q\\ \text{s.t.}& y\in Y\subseteq {\mathbb{R}}^n\\ & \text{Benders optimality cuts}\\ & \text{Benders feasibility cuts}\\ & q\text{ free}\end{array}$$
   其中，$q$是辅助变量，是为了衡量目前得到的$\text{Benders optimality cuts}$和$\text{Benders feasibility cuts}$为原始MIP的目标函数中关于$x$的部分，即$cx$的下界。在MIP的最优解中$cx$的取值$c{x}^{\ast }$，是等于在Benders分解中，得到最优解时候对应的$q^{\ast }$的。

这里，$\text{Benders optimality cuts}$的具体形式为
( α p i ) T ( b − B y ) ⩽ q , ∀ i ∈ P

$$\begin{aligned} (\alpha_p^i)^T(b-By) \leqslant q, \quad \forall i \in \mathcal{P} \end{aligned}$$ (αpi​)T(b−By)⩽q,∀i∈P​
其中，$\mathcal{P}$表示迭代过程中找到的子问题的对偶问题$\text{Dual SP}$的极点的集合，${\alpha }_p^i$就表示$\text{Dual SP}$的一个极点。

$\text{Benders feasibility cuts}$的具体形式为
( α r i ) T ( b − B y ) ⩽ 0 , ∀ i ∈ R

$$\begin{aligned} (\alpha_r^i)^T(b-By) \leqslant 0, \quad \forall i \in \mathcal{R} \end{aligned}$$ (αri​)T(b−By)⩽0,∀i∈R​
其中，$\mathcal{R}$表示迭代过程中找到的子问题的对偶问题$\text{Dual SP}$的极射线(extreme ray)的集合，${\alpha }_r^i$就表示$\text{Dual SP}$的一个极射线。

2. 求解主问题，得到$y$的值$\bar{y}$。并构建子问题。这里有2种等价的操作。

操作方法1：
根据$\bar{y}$构建子问题的对偶问题$\text{Dual SP}$
Dual SP: Q ( y ) = max ⁡ α α T ( b − B y ˉ ) s.t. A T α ⩽ c α  free

$$\begin{aligned} \text{Dual SP:} \quad Q(y) = \underset{\alpha}{\max} \quad & \alpha^T (b-B\bar{y}) \\ \text{s.t.} \quad &A^T\alpha \leqslant c \\ & \alpha \text{ free} \end{aligned}$$ Dual SP:Q(y)=αmax​s.t.​αT(b−Byˉ​)ATα⩽cα free​
优点：得到对偶变量方便，直接就是决策变量的取值。
缺点：需要写出子问题的对偶。

操作方法2：
这里，也可以不构建$\text{Dual SP}$，直接只构建$\text{SP}$即可。即，只需要构建：
SP: Q ( y ) = min ⁡ x c x s.t. A x = b − B y ˉ x ⩾ 0

$$\begin{aligned} \text{SP:} \quad Q(y) = \underset{x}{\min} \quad &cx \\ \text{s.t.} \quad &Ax = b - B\bar{y} \\ & x\geqslant 0 \end{aligned}$$ SP:Q(y)=xmin​s.t.​cxAx=b−Byˉ​x⩾0​
优点：不用写出子问题的对偶，建模方便。
缺点：CPLEX需要用IloCplex.getDuals()获取约束的对偶变量，需要用IloCplex.dualFarkas()函数获得约束的极射线； Gurobi中需要使用Constr.Pi获取约束的对偶变量，需要用Constr.FarkasDual获取约束的极射线。其实这个也不算啥缺点。

3. 求解子问题，获得对偶变量，构建Benders optimality cut和Benders feasibility cut。这里也有2种等价的操作。

操作方法1：
求解$\text{Dual SP}$，直接得到$\alpha$。

• 如果子问题的对偶$\text{Dual SP}$存在有界的可行解，则构建Benders optimality cut：
  $$\begin{array}{r}(\alpha {)}^T(b-By)⩽q\end{array}$$
• 如果子问题的对偶$\text{Dual SP}$无界，则构建Benders feasibility cut
  $$\begin{array}{r}(\alpha {)}^T(b-By)⩽0\end{array}$$

操作方法2：
直接求解$\text{SP}$的原始形式，然后用求解器得到$\text{SP}$的约束的对偶变量${\alpha }_p$或者极射线${\alpha }_r$。

• 如果子问题本身$\text{SP}$存在有界的可行解，则构建Benders optimality cut：
  $$\begin{array}{r}(\alpha {)}^T(b-By)⩽q\end{array}$$
• 如果子问题本身$\text{SP}$无可行解，则构建Benders feasibility cut
  $$\begin{array}{r}(\alpha {)}^T(b-By)⩽0\end{array}$$

4. 将步骤3构建好的Benders optimality cut和Benders feasibility cut添加到$\text{MP}$中。求解$\text{MP}$，并更新全局的上界$UB$和全局的下界$LB$。如果$UB=LB$，则算法停止，获得最优解，否则，循环2-4步。

这里需要说明，全局的$UB$和$LB$的计算方法如下：
U B = f y + Q ( y ) L B = f y + q

$$\begin{aligned} &UB = fy + Q(y) \\ &LB = fy + q \end{aligned}$$ ​UB=fy+Q(y)LB=fy+q​
这是为什么呢？因为$fy+Q(y)$，是给定了$y=\bar{y}$之后，求解$\text{SP}$，得到$x$的值$\bar{x}$，因此$(\bar{x},\bar{y})$是原问题MIP的一个可行解。$\min$问题的任意一个可行解，一定是原问题的$UB$，因此
$$\begin{array}{rl}& UB=fy+Q(y)\end{array}$$
。而$fy+q$是忽略了$x$的影响，只是添加了一部分的Benders optimality cut和Benders feasibility cut，并没有将所有的Cuts都加回来。因此它是原问题的一个松弛版本，因此可以为原问题提供一个下界。因此
$$\begin{array}{rl}& LB=fy+q\end{array}$$
当$UB=LB$时，我们就得到了全局最优解，此时正好满足
$$\begin{array}{rl}& fy+q=LB=UB=fy+Q(y)\end{array}$$
也就是
$$\begin{array}{r}q=Q(y)\end{array}$$

因此，比较简洁的判断Benders Decomposition是否得到最优解的办法，就是判断主问题中的$q$的取值，是否和子问题的目标函数$Q(y)$相同。

但是大家需要知道为什么是$q=Q(y)$。背后的原因，在上面的文字中已经给出了详细的理论解释。

直观上来讲，Benders Decomposition的总体框架如下。

基于此，我们给出Benders Decomposition的伪代码，方便后面的代码实现（注意，伪代码是基于$\min$问题的，如果是$\max$问题，可以相应做对应变化即可）：

Benders Decomposition的完整例子+代码实现

问题描述

文献[2]中给出了一个详细的例子，但是我觉得本文介绍的例子也许更适合来结合代码和可视化工具进行介绍。

该例子参考自网址：https://www.youtube.com/watch?app=desktop&v=vQzpydNOWDY

• Ray (Jian) Zhang. Learn Optimization Visually: Benders Decomposition. url: https://www.youtube.com/watch?app=desktop&v=vQzpydNOWDY

具体描述如下（有改动）：

小王拥有1000元现金流。他需要优化自己的理财计划，使得自己的收益最大化。他有以下两类选择。

1. 小王可以将钱存入银行储蓄账户，储蓄账户的固定收益率为4.5%。但是储蓄账户只允许存储整数数量的资金。
2. 小王还可以购买基金。假设一共有10个基金可选，分别为$F_1,F_2,\cdots ,F_{10}$，每个基金的收益率分别为1%, 2%, … 10%。每种基金最多购买100元。

问：小王应该如何投资，才能最大化自己的收益。

本问题的最优解为：

小王应该分别购买100元的$F_5,F$