优化 | Benders Decomposition
编者按： Benders算法虽然存在一定的弊端，但是Ta无可置疑的成为了运筹学史上最经典的算法之一。本篇文章从基础切入，直观的介绍了Benders算法的内核。
Benders Decomposition
针对由于规模过大或者结构过于复杂而无法直接整体求解的数学模型，可以将原问题分解成多个简单、能直接求解的小规模问题，再结合小规模问题求解得到的信息得到模型的精确解或近似精确解。Benders Decomposition就是这样的一种方法，它常用于解决随机规划问题和混合整数线性规划问题，是由Jacques F. Benders于1962年首先提出。
1. 基础知识回顾
1.1 线性规划的几何解释
对于一个线性规划问题，其 标准型 为
(LP) \quad Min \ \mathbf{c}^{T} \mathbf{x} \\ \quad s. t .\ \ A x=\mathbf{b} \\ \quad\quad x \geq 0 \\
可行域 为
P=\left\{\mathbf{x} \in \mathbf{R}^{n} \mid \mathbf{A} \mathbf{x}=\mathbf{b}, \mathbf{x} \geq \mathbf{0}\right\} \\
若 x 为该问题的可行解，则有 x \in P .
(1) 线性规划中的等式约束的几何意义
我们知道二维空间中， 2x+y=0 这样的一个等式约束表示一条直线，三维空间中一个等式约束表示一个平面，那么n维空间中等式约束表示什么呢？
引入**超平面（hyperplane）**定义：
对于一个向量 \mathbf{a} \in \mathbf{R}^{n} 且 \mathbf{a} \neq \mathbf{0} 和一个常数 \beta \in \mathbf{R} ，定义
H= \{x \in \mathbf{R}^{n}\left|\mathbf{a}^{T} \mathbf{x}=\beta \right \} \\
为一个超平面。
一个线性的等式约束即为一个超平面，它可以把整个空间分割成两个封闭半空间（halfspace）
\left\{\mathbf{x} \in \Re^{n} \mid \mathbf{a}^{T} \mathbf{x} \leq \beta \right\} \\ \left\{\mathbf{x} \in \Re^{n} \mid \mathbf{a}^{T} \mathbf{x} \geq \beta \right\} \\
(2) 线性规划可行解集合的几何意义
线性规划的可行域为凸集（convex set），即可行域内任意两点的连线都在可行域内，有限个线性约束所构成的可行域可以定义为一个 多面体（polyhedron） ，这个多面体是由 有限个封闭半空间的交集 所形成的。如果这个多面体是非空且有界的则可被称为一个 多胞形（polytope） 。
一个线性规划模型可能是不可行的（没有可行解）、无界的（有可行解却没有有限最优解），或是有有限最优解，对应的多面体分别为空、非空且无界和非空且有界。
我们可以用**极点（extreme point） 和 极方向（extreme direction or extreme ray，也叫极射线）**来描述多面体的特征，下面简单介绍一下极点和极方向。
- 极点对于求解优化问题是非常重要的，本科学习运筹学的时候，我们了解到如果线性规划问题是有最优解的，那么最优解肯定是在至少一个几何顶点（也是极点）中产生，而常用的单纯形法就是在极点上搜索最优解，极点对应 基本可行解（basic feasible solution） ，其实基本可行解就是极点的代数表达形式。
- 极方向的定义是如果一个向量 \mathbf{d}(\neq \mathbf{0}) \in \mathbf{R}^{n} ，满足
\left.\left \{\mathbf{x} \in \mathbf{R}^{n}\right| \mathbf{x}=\mathbf{x}^{0}+ \theta \mathbf{d} . \theta \geq 0\right\} \subset P, \ \ 对于所有\mathbf{x}^{0} \in P \\
那么向量 d 是多面体 P 的一个极方向，这里的 \theta 为步长。代入标准型的约束条件，若 d(\not=0) 是多面体的一个极方向，则有 Ad=0 且 d \geq 0 。对于一个无界的线性规划问题，可行域中的点沿着极方向移动，不管移动的步长有多大，仍在可行域内。若多面体是有界的，那么该多面体肯定不存在极方向，向量 d 为0；若多面体是无界的，那么该多面体存在极方向，向量 d 不为0。
如果我们可以找到线性规划所有的极点，那么加上极方向便可表示可行域中的所有点。可行域中的每个点都可以用其 极点 和 极方向 的线性组合表示出来（ Resolution theorem ），其中所有系数均为非负，且极点的系数之和为1，即
\mathbf{x}=\sum_{i=1}^{k} \lambda_{i} \mathbf{x}^{i}+ \theta \mathbf{d} \ \ 其中\lambda_{i} \geq 0, \theta \geq 0, \sum_{i=1}^{k} \lambda_{i}=1 \\
对于一个有可行解的最小化线性规划问题，其可行域为一个非空的多面体，它要么是有界的要么是无界的，不管是否有界，该问题的最优目标函数值要么至少在一个极点处取得，要么为负无穷，也就是该线性规划求解结果只有两种情况：一是多面体的极点，二是多面体的极方向。
根据 Resolution theorem ，将 \mathbf{x} 代入目标函数中，得到
\mathbf{c^{T}x}= \sum_{i=1}^{k} \lambda_{i} c^{T} \mathbf{x}^{i}+ \theta c^{T} \mathbf{d} \\
若多面体 P 存在极方向 d ，且满足 c^{T}d < 0 （说明向量 -c 和向量 d 之间的夹角小于90°），则 P 是无界的，且沿着 d 有 z \rightarrow-\infty ； 反之也是正确的，若最小化线性规划问题的目标函数 z \rightarrow-\infty ，那么多面体 P 一定存在满足 c^{T}d < 0 的极方向 d 。
1.2 线性规划的对偶理论
-
弱对偶理论Weak duality theorem
若 x 是原问题（最小化问题）的可行解， y 是相应的对偶问题的可行解，则有
\mathbf{y}^{T} \mathbf{b} \leq \mathbf{c}^{T} \mathbf{x} \\ 就是说，对于一个最小化问题，在有解的情况下，其对偶问题可行解的目标函数值可为该问题的目标函数值提供 下界 ；而对于一个最大化问题，其对偶问题可行解的目标函数值可为该问题的目标函数值提供 上界 。
我们知道线性规划问题可能是无解的或者无界的，弱对偶理论为我们分析无界和无解的情况提供工具。(1) 假设一个最小化问题 无界 ，弱对偶定理告诉我们对偶可行解的目标函数值可为原问题提供下界，而此时原问题的解无界，说明没有这样的下界存在，也就是说没有对偶可行解存在，那么 对偶问题无解 。(2) 假设一个最小化问题 无解 ，我们把原问题看成是其对偶问题（最大化问题）的对偶，从弱对偶定理我们可以知道原问题无法为对偶问题提供上界，即无法限定对偶问题的可行解范围，那么 对偶问题是无界或是无解的 。弱对偶定理说明了原始和对偶可行解的目标函数值互相为对方的界。
-
强对偶理论Strong duality theorem
如果一个线性规划问题或者其对偶问题有有限最优解，那么该线性规划及其对偶问题均有有限最优解，且两者最优目标函数值相等。
2. Benders Decomposition原理浅析
Benders Decomposition的基本思想是通过将求解难度较大的原问题分解为两个更容易求解的问题master problem和subproblem，类似于列生成（Column Generation）方法，不同的是列生成方法通过求解subproblem产生column添加至master problem，而Benders Decomposition通过求解subproblem产生cut添加至主问题（类似row/ constraint generation）。如何产生什么样的cut是我们关心的问题。
2.1 问题分解及cut的产生
Benders Decomposition多适用于求解如下形式的混合整数规划（mixed integer programming，下文简写为MIP）模型：
(MIP) \quad \min \quad c^{T} x+f^{T} y \\ s. t. \quad \mathbf{A x}+\mathbf{B y} = \mathbf{b}\\ \mathbf{x} \geqslant \mathbf{0}, \quad \mathbf{y} \geqslant 0 \quad 且为整数 \\
其中 x 为连续变量， y 为整数变量。由于原问题模型中决策变量同时包含连续变量和整数变量，并且不同类型的变量通过约束条件 \mathbf{A x}+\mathbf{B y} =\mathbf{b} 联系在一起，增加了问题的求解难度。我们要想办法将离散变量与连续变量分离。
固定 y 的值为非负整数 \hat{y} ，则对应一个关于连续变量 x 的线性规划问题，看成是一个关于 \hat{y} 的函数 g(\hat{y}) ，则有：
(M1) \quad g(\hat{y}) =\min \ c^{T} x \\ A x =b-B \hat{y} \\ x \geq 0 \\
那么原MIP则变为如下形式：
(M2)\quad \min \quad f^{T} y + g(\hat{y}) \\ s. t. \quad \quad \mathbf{y} \geq 0 \quad 且为整数 \\
定义(M1)为primal subproblem，定义(M2)为master problem。
这样我们便将求解困难的MIP分解成相对容易求解的 master problem 和 subproblem。
通过观察我们发现subproblem的右端项 b-B \hat{y} 是随 y 的不同取值而变化的，这样每次求解 primal subproblem 的线性规划，其可行域是变化的，而它的对偶问题，可行域是不变的，那么每次迭代只需更新目标函数的价格系数 b-B \hat{y} 。写出 primal subproblem 的对偶问题 dual subproblem(M3)，假设 subproblem 的约束 A x =b-B \hat{y} 的对偶变量为 \pi ，则有：
(M3) \quad \max \pi^{T}(\boldsymbol{b}-\boldsymbol{B} \boldsymbol{\hat{y}}) \\ s.t. \quad \boldsymbol{\pi}^{T} \boldsymbol{A} \leq \boldsymbol{c} \\ \\
subproblem 的目标是找到确定的 \hat{y} 下最优的 x （或者 \pi ），可以看到 dual subproblem 的多面体 P 不随变量 y 的取值变化而变化，它的求解情况有三种可能：无解、有可行解但无界（无界解）、有可行解且有界（有界解）。
- 当 dual subproblem 无解时，多面体 P 为空集，由弱对偶定理可知primal subproblem无解或无界，那么master problem也无解或无界，原问题MIP同样也是无解或无界。
- 当 dual subproblem 的解为有界解时，多面体 P 不为空，由弱对偶定理可知 primal subproblem 有有界解。
- 当 dual subproblem 存在无界解时，多面体 P 不为空，由弱对偶定理可知 primal subproblem 无解。
Benders Decomposition就是依据 dual subproblem 求解所得到的解的情况来产生添加到master problem的cut，每添加一个cut 都会将原来的多面体割掉一部分，具体分析如下：
-
当 dual subproblem 存在无界解时， primal subproblem 无解，说明给定的
\hat{y}
不合理，需要添加cut割掉
\hat{y}
，最简单的方法就是找到
\hat{y}
所违反的约束，加入master problem。
由上文提到的 Resolution theorem 可知，对于最大化（最小化）线性规划问题存在无界解，那么多面体一定存在极方向 d 使得目标函数趋近于无穷大（小），并且极方向 d 满足条件 c^{T}d > 0 （ c^{T}d < 0 ）。
假设多面体 P 的极方向集合为 \Omega_{d} ，那么存在 \boldsymbol{w} \in \Omega_{d} ，使得 \boldsymbol{\pi}_{w}^{T}(\boldsymbol{b}-\boldsymbol{B} \boldsymbol{\hat{y}})>0 ，向 master problem 中添加cut：$\boldsymbol{\pi} {w}^{T}(\boldsymbol{b}-\boldsymbol{B} \boldsymbol{\hat{y}}) \leq0, \forall \boldsymbol{w} \in \Omega {d} 便割掉了 \hat{y}$ ，添加的约束就是 feasibility cut 。
-
假设多面体
P
的极点集合为
\Omega_{p}
，当 dual subproblem 的解为有界解时，最优解在至少一个极点处
{\pi}_{w}
取得，其中
\boldsymbol{w} \in \Omega_{p}
，primal subproblem 的解也为有界解。根据
Resolution theorem
，多面体
P
不存在使得
\boldsymbol{\pi}_{w}^{T}(\boldsymbol{b}-\boldsymbol{B} \boldsymbol{\hat{y}})>0
的
{\pi}_{w}
，其中
\boldsymbol{w} \in \Omega_{d}
，即
\boldsymbol{\pi}_{w}^{T}(\boldsymbol{b}-\boldsymbol{B} \boldsymbol{\hat{y}}) \leq0
。根据强对偶理论，这种情况下dual subproblem 和 primal subproblem 的最优目标函数值相等，假设为
\eta
，可以通过向 master problem 添加cut：
\boldsymbol{\pi}_{w}^{T}(\boldsymbol{b}-\boldsymbol{B} \boldsymbol{\hat{y}}) \leq \eta,\forall \boldsymbol{w} \in \Omega_{p}
来继续优化master problem，改进 master problem 的目标函数值，朝更小的目标前进，添加的约束就是
optimality cut
。
- 总结如下表
2.2 Benders Decomposition算法过程
Step 1 . 初始化：设置MIP下界为 -\infty ，上界为 +\infty ，固定整数变量 y= \hat{y} \in Y ，代入MIP约束 \mathbf{A x}+\mathbf{B y} = \mathbf{b} ，则 f^{T} \hat{y}、 \mathbf{B \hat{y}} 为常数，将MIP分解为 relaxed master problem 和 primal subproblem，由于master problem没有cuts（相当于松弛掉了很多约束后的一个模型，relaxed master problem 的最优目标函数值为MIP提供下界）；
Step 2 . 求解 dual subproblem，最优解为 \pi^{*}
- 若 \pi^{*} 为无界解，向 relaxed master problem 添加 feasibility cut
- 若 \pi^{*} 为有界解，向 relaxed master problem 添加 optimality cut
Step 3 . 求解 relaxed master problem，更新MIP上下界并得到最优解 y^{*} ，加入 dual subproblem，回到Step 2 或进入Step4
- 若添加 feasibility cut，求解 relaxed master problem最优目标函数值为MIP下界
- 若添加 optimality cut，求解 relaxed master problem最优目标函数值为MIP上界
Step 4. 若满足终止条件（比如上下界之间的gap小于0.01），则停止迭代。
Benders Decomposition的伪代码：
\begin{align*} & Initialize \ \ y^{*}, U B=+\infty, L B=-\infty \\ & While \ \ U B-L B>\varepsilon \\ & Solve \ the \ Dual\ Subproblem, \pi^{*} \\ & If \ \ \pi^{*} \ is\ unbounded \\ & Add\ cut \ [b-B y]^{T} \pi^{*} \leq 0 \\ & Else \ \ L B:=\max \left\{L B, f^{T} y^{*}+\left[b-B y^{*}\right]^{T} \pi^{*}\right\} \\ & Add \ cut \ z \geq f^{T} y+[b-B y]^{T} \pi^{*} \\ & End\ If \\ & Solve \ Master\ problem \ \max _{y}\{z \mid c u t s, y \in Y\} \\ & U B:=z^{*} \\ & End \ While \end{align*} \\
简单的 master problem 和 subproblem可以直接用Gurobi、Cplex等求解器求解，对于复杂的 master problem 和 subproblem 用求解器直接求解往往很困难，需要采用比如分支定界、分支切割等优化求解方法来求解。
参考文献
- [1] Bertsimas D, Tsitsiklis J N. Introduction to linear optimization[M]. Belmont, MA: Athena Scientific, 1997.
- [2] Rardin R L, Rardin R L. Optimization in operations research[M]. Upper Saddle River, NJ: Prentice Hall, 1998.
原创声明
作者：莫斑炜
责任编辑：苏向阳
审核编辑：阿春
微信编辑：玖蓁
知乎编辑：乐陶
本文由『运筹OR帷幄』 @运筹OR帷幄 原创发布，转载请至同名公众号后台获取转载须知
