最近在学习Benders Decomposition，看了不少的文章博客，但我感觉写的都不是很清楚，细节地方解释的也不够到位，于是就去读了英文论文原文，发现原文写的还是相当清楚的。本文也基于那篇论文，对Benders Decomposition进行总结。

Benders Decomposition由Jacques F. Benders提出，用来求解混合整数规划问题(MIP:Mixed Integer Programming problem)，即同时包含整数和连续变量的极值问题，换句话说就是线性规划和整数规划结合起来的数学规划问题。

原理介绍：Benders分解算法的精妙之处在于引入了复杂变量（complicating variables），当这些变量固定后，剩下的优化问题（通常称为子问题）变得相对容易。在MIP中，先把复杂变量（整数变量）的值固定，则问题成为一般的线性规划问题，当然，这个线性规划问题是以复杂变量为参数的。在Benders设计的算法里，利用割平面的方式将主问题（以子问题的解为参变量）的极值和使子问题（线性规划问题）有可行解的参变量值的集合很恰当地表达了出来。

首先，给出IMP原问题的矩阵表述形式：（1）

 Minimize  c T x + f T y  subject to:  A x + B y = b x ≥ 0 , y ∈ Y ⊆ R q \begin{array}{l} \text { Minimize } c^{T} x+f^{T} y \\ \text { subject to: } A x+B y=b \\ \qquad x \geq 0, y \in Y \subseteq \mathbb{R}^{q} \end{array}  Minimize cTx+fTy subject to: Ax+By=bx≥0,y∈Y⊆Rq​

这里的 x x x和 y y y分别是 p p p和 q q q维向量（列）， Y Y Y是 y y y的可行域， A A A、 B B B是对应维度的矩阵， c c c、 b b b、 f f f也是对应维度的向量。

这里，我们把 y y y当成复杂变量（ y y y是整数好像并不复杂，在此我的理解是 y y y的约束更细致，所以“复杂”），当 y y y值固定时，原问题就变成了普通的线性规划问题。基于此，我们可以把问题进行分解：

主问题(master problem):（2）

 Minimize  f T y + q ( y )  subject to:  y ∈ Y \begin{array}{l} \text { Minimize } f^{T} y+q(y) \\ \text { subject to: } y \in Y \end{array}  Minimize fTy+q(y) subject to: y∈Y​

子问题(subproblem):（3）

 Minimize  c T x  subject to:  A x = b − B y x ≥ 0 \begin{array}{l} \text { Minimize } c^{T} x \\ \text { subject to: } A x=b-B y \\ \qquad x \geq 0 \end{array}  Minimize cTx subject to: Ax=b−Byx≥0​

这里， q ( y ) q(y) q(y)是子问题的最优目标函数值( the optimal objective function value)。对于任意给定的 y ∈ Y y\in Y y∈Y，子问题是一个线性规划问题。可以得到，如果子问题无界（unbounded），那么主问题也必定无界，则原问题无最优解；在子问题有界( boundedness)的情况下，我们可以通过求解子问题的对偶问题来计算 q ( y ) q(y) q(y)。

所以，我们对子问题引入对偶变量 α \alpha α，那么子问题的对偶问题就是：（4）

M a x i m i z e   α T ( b − B y ) s u b j e c t t o : A T α ≤ c α   u n r e s t r i c t e d . \begin{array}{l} Maximize \space \alpha^{T}(b-B y)\\ subject to: A^{T} \alpha \leq c \\ \alpha \space unrestricted.\\ \end{array} Maximize αT(b−By)subjectto:ATα≤cα unrestricted.​

可以发现，对偶问题的可行域（feasible region）并不依赖于 y y y， y y y只影响目标函数值。因此，对于给定的 y y y，如果对偶问题可行域为空，那么原问题无界或可行域为空（此时任意 y ∈ Y y\in Y y∈Y原问题都为空）。我们假设对偶问题可行域不为空，我们可以枚举可行域所有的极点（extreme points）： ( α p 1 , . . . . , α p I ) (\alpha_p^1,....,\alpha_p^I) (αp1​,....,αpI​)，极射线(extreme rays)： ( α r 1 , . . . , α r J ) (\alpha_r^1,...,\alpha_r^J) (αr1​,...,αrJ​)， I I I和 J J J分别指的是极点和极射线的数量。然后，对于一个给定的 y ^ \hat{y} y^​向量，要求解对偶问题，可以通过检测：（i） ( α r j ) T ( b − B y ^ ) > 0 (\alpha_r^j)^T(b-B\hat{y})>0 (αrj​)T(b−By^​)>0对于所有极射线 α r j \alpha_r^j αrj​是否成立，如果成立则对偶问题无界且原问题无解；（ii）是否存在一个极点使对偶问题的目标函数值 ( α p i ) T ( b − B y ^ ) (\alpha_p^i)^T(b-B\hat{y}) (αpi​)T(b−By^​)最大，若存在，则原问题和对偶问题都有有限的最优解（finite optimal solutions）。

基于上述理论，对偶问题可以重写成以下形式：（5） KaTeX parse error: Double superscript at position 140: …\alpha_{p}^{i} ^̲{T}(b-B y) \leq…

从上式可以发现，只有 q q q一个变量，但是存在大量约束。我们可以用 q ( y ) q(y) q(y)来替代主问题，这样就可以得到只有 q q q和 y y y变量的形式：（6） M i n i m i z e f T y + q s u b j e c t   t o : ( α r j ) T ( b − B y ) ≤ 0 ∀ j = 1 , … , J ( b ) ( α p i ) T ( b − B y ) ≤ q ∀ i = 1 , … , I ( c ) y ∈ Y , q   u n r e s t r i c t e d . Minimize f^{T} y+q\\ subject \space to: \\ \left(\alpha_{r}^{j}\right)^{T}(b-B y) \leq 0 \quad \forall j=1, \ldots, J(b)\\ \left(\alpha_{p}^{i}\right)^{T}(b-B y) \leq q \quad \forall i=1, \ldots, I (c)\\ y \in Y, q\space unrestricted. MinimizefTy+qsubject to:(αrj​)T(b−By)≤0∀j=1,…,J(b)(αpi​)T(b−By)≤q∀i=1,…,I(c)y∈Y,q unrestricted.

由于对偶表达式（ dual formulation）存在大量极射线和极点，要生成上述所有的约束是不现实的。所以Benders分解算法使用约束条件的子集，即求解松弛主问题（relaxed master problem）。开始，初始松弛主问题中无约束，在Benders算法求解过程中不断向松弛主问题中加入式（6）中约束条件的某一个，即加入有效的切平面（cut）。通过求解松弛主问题，我们可以得到一个候选最优解 ( y ∗ , q ∗ ) (y^*,q^*) (y∗,q∗)，然后将 y ∗ y^* y∗代入对偶子问题（4）中求解计算 q ( y ∗ ) q(y^*) q(y∗)值，如果子问题的最优解 q ( y ∗ ) = q ∗ q(y^*)=q^* q(y∗)=q∗，则算法停止。如果对偶问题无解，则在松弛主问题中可以加入（6b）类型的约束，然后求解新的松弛主问题。（6b）类型的约束称之为Benders feasibility cuts。如果对偶子问题的最优解 q ( y ∗ ) > q ∗ q(y^*)>q^* q(y∗)>q∗，则在松弛主问题中可以引入（6c）类型的约束，然后求解新的松弛主问题。（6c）类型的约束称之为Benders optimality cuts。在每次迭代过程中都可以生产某一类型的约束，由于 I I I和 J J J是有限的，故可以保证在有限次迭代过程后得到最优解。

初始化： { y = 有效的整数解 LB = -inf UP = inf } while(UB - LB > 0) { 求解子问题 m i n u { f T y ^ + ( b − B y ^ ) T α ∣ A T α ≤ c , α ≥ 0 } min_u\{f^T\hat{y}+(b-B\hat{y})^T\alpha|A^T\alpha\leq c,\alpha\geq0 \} minu​{fTy^​+(b−By^​)Tα∣ATα≤c,α≥0}

​ if 无界 {

​ 得到无界的极射线 α r \alpha_r αr​

​ 把割 ( b − B y ) T α r ≤ 0 (b-By)^T\alpha_r\leq 0 (b−By)Tαr​≤0加入到主问题 } ​ else { ​ 得到极点 α p \alpha_p αp​

​ 把割 z ≥ f T y + ( b − B y ) T α p z\geq f^Ty+(b-By)^T\alpha_p z≥fTy+(b−By)Tαp​加入到主问题

​ U B : = min ⁡ { U B , f T y ^ + ( b − B y ^ ) T α p } U B:=\min \left\{U B, f^{T} \hat{y}+(b-B \hat{y})^{T} \alpha_p\right\} UB:=min{UB,fTy^​+(b−By^​)Tαp​} } ​ 求解主问题

​ min ⁡ y { z ∣  cuts,  y ∈ Y } \min _{y}\{z \mid \text { cuts, } y \in Y\} miny​{z∣ cuts, y∈Y}

​ LB = z

}

本来还想最后加个源代码实战教学一下，但是感觉这个代码我一时半会儿还写不出来，以后写完再放出来吧。

把东西写出来还是相当麻烦耗费时间的，有什么问题可以在评论中指出来交流一下！