Benders分解

参考视频：068-线性规划（十三）：Benders分解-哔哩哔哩

1. Benders分解简介

1.1 什么是Benders分解

Benders分解算法是由提出它的Jacobus Franciscus (Jacques) Benders教授的名字命名的，核心思想是将两类变量（一类为整数变量，另一类为连续变量）分开求解，然后再合起来。

举个例子，这就像小鸡炖蘑菇：

传统的做法是把小鸡和蘑菇一锅炖出来。事实上，小鸡比较难炖，蘑菇比较好炖。考虑到这一点，Benders分解的思想就是把小鸡（复杂变量）和蘑菇（简单变量）分开炖，但这两锅又要有联系：因为蘑菇要用鸡汤来炖，在炖的过程中我们又要尝尝蘑菇鸡汤的咸淡，再决定要不要往鸡汤里放更多的盐。也就是说，炖鸡的那口锅要给炖蘑菇的那口锅提供鸡汤，而炖蘑菇的那口锅要给炖鸡的那口锅反馈咸淡。重复上面的过程，直到分开炖的口味和一锅炖出来的口味相差无几，我们就把两口锅的东西合起来，就得到了跟传统做法一样原汁原味的小鸡炖蘑菇(๑´ڡ`๑)

1.2 学习Benders分解的前置知识

• Cutting plane（割平面）

• Vaild inequalities（有效不等式）：也就是割平面的过程中，利用有效不等式，割去一部分内容物，但没有割去最优解（保持最优特性）

• Row generation（行生成）：行生成的概念

• Duality theory（对偶理论）：最核心的知识，详情请见之前的博客大规模线性规划的对偶 - 码头牛牛 - 博客园

首先来介绍一下Benders分解的初步逻辑：

• 如图，Benders分解会把一个优化问题分解为一个主问题(Master Problem)和一个子问题(Sub-Problem)分别进行求解。

  • Master Problem (MP)：只剩下复杂变量（在图中指的是整数变量\(y\)）和用来传递信息的变量（在图中指的是变量\(q\)）的优化问题

  • Sub-Problem (SP)：包含简单变量的优化问题（图中指的是连续变量\(x\)）

• 在求解过程，MP会先指定一个可行的\(y\)（相当于控制变量），传给SP。

• SP会反馈给MP这个\(y\)是否合适。如果合适的话会利用\(y\)的信息帮助找到最优解；如果\(y\)提供的信息无效或者说是不可行的，那就要避免MP下次还提供这样的信息，也就是把这样的信息cut掉。也就是：

然后再简单介绍一下行生成（有关列生成内容在这里：列生成算法 - 码头牛牛 - 博客园）：

• 行生成关注的是约束，因为约束是以行的形式存在的；列生成关注的是基，因为基是以列的形式存在的。

• 行生成是从一个隐藏的空间里，不断地去拿有效的约束（有效不等式），并把它们放到松弛的MP里；列生成是从对偶子问题里面，把一个reduce cost比较好的基反馈给MP。

• 比较典型的行生成方法有：割平面、User Cut等等，下面介绍其中两个：

  • User cut：就是一个提前加的割，作用是减小搜寻范围。在我们对问题有初步把握的基础上，我们可以人工加入几个割平面 (User cut)，来割掉一些多余的搜寻空间，以加速搜索进程

  • Lazy cut：当需要的时候就添加，不需要的时候就不添加。作用是使解变得更加可行。它的思想和Benders分解中加cut的思想是一致的。

  • 如下面第一张图，黄色部分是松弛后的可行域。当我们比较了解这个问题的时候，就可以事先加入User Cut，来缩小这个黄色的搜寻范围；如下面第二张图，有些问题为了方便求解，需要把一些约束给拿掉，拿掉约束后，它会添加一些不可行的点进来，也就是图中灰色的点。这导致了求解的过程中可能会碰到这些点，如果碰到了发现不可行，我们就会添加一个Lazy Cut进来，把这个点砍掉。所以说Lazy cut就是使用到的时候才会拿出来用，而Benders分解的思想就是不断地向松弛主问题添加Lazy Cut。

2. Decomposition process

如图，给出一个具有两类变量的问题，\(x\)是简单的连续变量，\(y\)是复杂变量，\(y\)的可行域是\(Y\)。然后我们把它分解成两类问题（图中红色的\(y\)是一个变量，绿色的\(y\)是一个给定的值）：

• SP：只包含\(x\)，在子问题中\(y\)是需要提前给定的。图中SP的目标函数中，\(f^Ty\)其实加不加都行，因为这是个给定常数。

• MP：包含变量\(y\)的主问题，它只含有复杂变量。

求解MP，为SP提供一个确定的\(y\)；SP会得到当前\(y\)的求解情况，并以割平面(Cutting plane) 的方式反馈给MP。

• 重点在于如何找到这个割平面，以指导MP下次这个\(y\)的生成。

• 要想找这个割平面，首先就要知道MP生成的\(y\)是如何影响SP的。只有知道输入的\(y\)如何影响SP，才能知道该反馈什么样的信息给MP，以指导MP生成更好的\(y\)。

• 所以在这里，割平面相当于一个路径，这个路径让我们知道什么会影响\(y\)的取值。

2.1 SP和MP之间的关系

传入的\(y\)如何影响SP? 为了研究这个问题，我们把MP和SP合并起来写：

我们对SP提供一个\(y\)，会出现以下两种情况：

• SP feasible，这种情况又可以分为三类：

  • Unique solution(唯一的可行解)

  • Infinite solutions(无穷多个可行解)：它是一个有界解，只不过可行解是一个平面或一个空间

  • Unbound solution(无界解)：也就是可行解是开放的，这就导致了目标函数值会变成一个无穷

• SP infeasible，说明\(y\)输入SP后，导致SP无可行域。此时这个节点就相当于一个叶子节点，搜索也就到此打住了。

针对SP得到的四种不同的情况（3种不同的解及不可行），MP会出现以下3种情况：

• Feasible Bounded：当SP存在Unique solution或Infinite solutions，MP就是一个可行有界的情况。

• Feasible unbounded：如果SP的解是Unbound solution，那么SP会朝着一个方向不断优化，一直到无穷无尽，这也会导致MP变得无穷无尽。

• Infeasible：如果SP infeasible，那MP也是infeasible

由此，我们也就建立了SP和MP之间的关系，也把它们之间的相互影响路径找出来了。

2.2 SP如何向MP传递信息（D-SP、SP、MP之间的相互影响路径）

2.1中SP和MP的相互影响路径有一个缺点：MP每次向SP传递一个\(y\)，SP都要重新计算它的可行域，所以这个路径我们一般不常用。

我们发现SP是一个线性规划问题，这时候我们尝试看看能不能用它的对偶问题(D-SP) 来建立这个路径：

• 其中，\(\alpha\)是对偶变量。

2.1已明确了SP和MP之间的相互影响路径，而如果想明确D-SP和MP之间的相互影响路径，我们可以先探究D-SP是如何影响SP的。

MP会传递一个\(y\)给SP，而我们发现，根据SP写出的D-SP中，这个\(y\)只存在于目标函数里，并没有出现在它的可行域（约束条件）。这是一个很好的性质：因为对于输入的\(y\)，我们不需要重新计算D-SP的可行域，也就意味着从迭代开始到最后，我们可以一直使用同一个可行域以及它的极点和极方向。

将\(y\)输入到D-SP中，它的可行域是固定的，它也会出现以下两种情况：

• Feasible (Unique solution / Infinite solutions / Unbounded solution)

• Infeasible

根据弱对偶性的推论，我们可以得到上述情况下SP的解会出现的情况：

• 如果D-SP的解是有界可行的，那么SP的解也是有界可行的。

• 如果D-SP的解是无界解，那么SP无可行解。

• 如果D-SP的解不可行，那么SP的解是无界解或不可行解。

然后根据2.1中SP和MP的关系，我们得到了D-SP、SP、MP之间的相互影响路径：

• MP要向D-SP提供一个\(y\)，在固定的\(y\)下求解D-SP，根据D-SP的解的情况，D-SP会反馈给MP这个\(y\)的好坏。

• 由于D-SP的可行域不受\(y\)的影响，所以如果D-SP无可行解的话，MP也没办法通过改变\(y\)来调节。也就是说，如果D-SP是Infeasible的话，那不管\(y\)是多少，它永远都是Infeasible，那这个相互影响路径也就没用了，这种情况是没有优化意义的（因为优化的结果要么是无解，要么解是一个无穷）。

根据图中的相互影响路径，我们下面分别展开讨论D-SP不同情况。

2.2.1 D-SP是Feasible Bounded的情况

如果D-SP的解是Feasible Bounded的，那么SP、MP的解也是Feasible Bounded的。

• 补充一下极点和极方向。比如\(P=\{(x_1,x_2)|x_1\geq0,x_2\geq0,x_1+x_2\geq1\}\).在这个例子中：

  • 极点 (extreme points)：\((1,0),(0,1)\)

  • 极方向(extreme ray)：\(d_1=(1,0),d_2=(0,1)\)

  • 我觉得极点比较好理解，极方向\(d_1、d_2\)都是一个向量，可以理解成从任意一个可行的\(x\)出发，它沿着向量\(d_1、d_2\)所表示的方向一直走，仍然留在可行域中。也就是说，对于任意可行解\(x\)，\(x+\lambda{d_1}\)或\(x+\lambda{d}_2\)（\(\lambda\geq0\)）都在可行域里。

根据对偶定理，在MP传入的这个\(\overline{y}\)下，SP的最小值=D-SP的最大值=SP、D-SP的最优值

换句话来说，SP的每个可行解都是D-SP的上界，我们只要找到D-SP最小的上界，就可以得到\(\overline{y}\)下SP、D-SP的最优解。那么有：

\[q(\overline{y})=min\quad{q} \]

\[\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad{s}.t.\quad(b-B\overline{y})^T\alpha^p_i\leq{q}\quad{i\in{I}} \]

• \(q\)：D-SP的一个上界。这里的目标函数是找一个最小的D-SP的上界（最小的\(q\)相当于D-SP目标函数的最大值）。

  • 需要注意的是，\(q\)相当于由SP的可行解算出来的目标函数值，而不是一个随便取的值。

  • 或者可以理解为，\(q\)大于或等于每一个极点\(\alpha^p_i\)算出来的D-SP目标函数值。在这里我个人的理解是，最后算出来加入MP的Optimal cut，它相当于\(max\{(b-B\overline{y})\alpha^p_i|i\in{I}\}\leq{q}\)

• \(\alpha^p_i\)：D-SP可行域的一个极点（也是一个可行解）。这里的约束条件相当于D-SP在某个极点处的目标函数值（这里为了方便直接把\(f^T\overline{y}\)去掉了，因为这是一个常数，并不会影响优化结果），小于等于它的一个上界。

• \(q(\overline{y})\)：就是在上图可行域的若干个极点中，挑一个最好的。

由此，我们把D-SP这个最大化问题形式，替换成了上面这个最小化问题的形式。我们把\(\quad(b-B\overline{y})^T\alpha^p_i\leq{q}\quad{i\in{I}}\) 称为Optimal cut.

---

下面是关于我对Optimal cut的理解。

Optimal cut的作用：把当前输入的\(y\)下SP的真实最优值信息传回 MP，防止 MP 把 \(q(y)\) 估得过小。

举一个例子便于理解Optimal cut的作用：

[原问题]

\[\min_{x,y}\quad{3y+x} \]

\[s.t.\quad{x+2y}\geq5 \]

\[x\geq0,y\in\{0,1,2,3\} \]

• 这个原问题的最优解是：\(y^*=0；x^*=5\)，最优值为\(5\)

[MP]

\[\min_{y\in\{0,1,2,3\}}\quad{3y}+q(y) \]

[SP]

\[q(y)=\min_x\{x:x\geq5-2y,x\geq0\} \]

• 在求解过程中，MP会基于\(y\in\{0,1,2,3\}\)这个范围，不断输入\(y\)到SP中。也就是说，对于SP，固定不同的\(y\)后，会有以下情况（下面情况都属于Feasible Bounded）：

  • \(y=0\)，解得：\(x^*=5\)，即\(q(0)=5\)

  • \(y=1\)，解得：\(x^*=3\)，即\(q(1)=3\)

  • \(y=2\)，解得：\(x^*=1\)，即\(q(2)=1\)

  • \(y=3\)，解得：\(x^*=0\)，即\(q(3)=0\)

我们现在先忽略SP，假设随便拿一个变量\(q\)来替代掉SP目标函数值的变量\(q(y)\)，那么，MP就会变为：

[不考虑SP的MP]

\[\min_{y\in\{0,1,2,3\}}\quad{3y}+q \]

这个优化问题目测一下就可以知道，肯定是 \(y\)和\(q\)越小越好，\(y\)最小是\(0\)，\(q\)最小可以取到\(-∞\)，但显然，原问题的最优解不是\(-∞\)

但实际上，受SP中\(x\)可行域的约束，\(y=0\)的时候，真实的\(q(y)=q(0)=5\)。也就是说，如果不考虑SP中那些关于\(x\)的约束，MP会倾向于把 \(q(y)\)（也就是上面的\(q\)） 压得尽可能小。

总结一下，MP实际上是不知道\(q(y)\)的取值是多少的，但因为它的目标函数是\(min\)，所以它会把\(q(y)\)的值尽可能压小。所以Optimal cut的作用是：告诉MP，当\(y\)取某个值时，\(q(y)\)（或者\(q\)）至少得等于这个SP的真实最优值。 然后通过上面Optimal cut的推导，使用对偶理论，我们又恰好得到了在输入的\(y\)下，\(q(y)\)的下界。

2.2.2 D-SP是Feasible Unbounded的情况

由于D-SP的可行域是固定不变的，也就是说，不同的\(y\)输进来，相当于在同一个可行域中换不同的的目标方向去做最大化。因此，若某个\(y\)对应的目标方向使目标函数能够沿可行域中的某个极方向无限增大，则D-SP就会出现 feasible unbounded 的情况，所以需要通过 feasible cut 将其切除。

MP传入一个\(\overline{y}\)到SP，这个\(\overline{y}\)导致了SP不可行和D-SP无界。feasible cut的作用就是把导致SP不可行和D-SP无界的\(\overline{y}\)从MP的可行域中去掉。

在D-SP是feasible unbounded的情况如下图：

• 此时，D-SP可行域里至少有一个可行解\(\alpha^0\)，并且存在某个方向\(\alpha^r_j≠0\)，使得对任意的\(\lambda\geq0\)，\(\alpha^0+\lambda\alpha^r_j\)仍然可行。随着\(\lambda\)的不断增大，目标函数\((b-B\overline{y})^T(\alpha^0+\lambda\alpha^r_j)\)不断增大，最后会变成一个无穷（这里目标函数是\(max\)，所以最后会变成正无穷）。

• \((b-B\overline{y})^T(\alpha^0+\lambda\alpha^r_j)=(b-B\overline{y})^T\alpha^0+\lambda(b-B\overline{y})^T\alpha^r_j\)

  • 其中，\(\alpha^0\)和\(\alpha^r_j\)都是一个已知的可行解和极方向，那么这个目标函数也就变成了一个关于 \(λ\) 的一元函数。

  • 对目标函数求导，可以得到它沿极方向\(\alpha^r_j\)走的增长率：\((b-B\overline{y})^T\alpha^r_j\)

    • 如果\((b-B\overline{y})^T\alpha^r_j>0\)，那么沿这个方向走，目标函数会越走越大，D-SP 无界，进而 SP 不可行。

    • 所以为了避免SP不可行，必须强制\((b-B\overline{y})^T\alpha^r_j\leq0\quad{j}\in{J}\)，这就是 feasible cut。

2.2.3 Optimal cut和Feasible cut的作用总结

Optimal cut：

• 用于D-SP可行有界的情况。Benders分解就是把原问题拆成了MP和SP，MP负责在\(y∈Y\)这个区域内生成\(y\)，在SP是Feasible Bounded的情况下，如果不考虑SP中\(Ax=b-By\)的约束，MP目标函数中的\(q(y)\)肯定是倾向于越小越好的。所以，在SP可行有界时，Optimal cut恰恰是考虑了SP中\(x\)的约束，限制了\(q(y)\)的下限。它的本质是为了让 MP 对目标值的刻画越来越准确。

Feasible cut：

• 用于D-SP可行但无界的情况。此时说明当前\(y\)会导致D-SP的目标值沿某个极方向不断增大到 \(+∞\)，进一步说明原来的SP不可行。因此，feasible cut 的作用是：把这类会导致SP不可行的\(y\)从 MP 的可行域中排除掉，避免MP以后再给出类似的无效解。

2.3 Benders 分解流程

2.3.1 松弛主问题

我们来回顾一下MP:

[MP]

\[\min_y\quad{q}(y)+f^Ty \]

\[s.t.\quad{y}\in{Y} \]

根据上面的Optimal cut和Feasible cut，对于MP传入的\(\overline{y}\)，D-SP也可以写成：

[D-SP(\(\overline{y}\))]

\[min\quad{q}+f^T\overline{y} \]

\[s.t.\quad(b-B\overline{y})^T\alpha^p_i\leq{q}\quad{i\in{I}} \]

或

\[min\quad{q}+f^T\overline{y} \]

\[s.t.\quad(b-B\overline{y})^T\alpha^r_j\leq{0}\quad{j\in{J}} \]

假设\(y\in\{y_1,y_2,y_3\}\)，那么结合D-SP，完整的MP可以写成（也就是我们不再把SP套入MP中，而是直接把每一种情况枚举到MP中）：

在现实应用中，\(y\)的取值是非常多的，这也导致了完整的MP的约束会很多（行很多），在这种情况下，我们用行生成的思想，先把上面那些cut全拿掉，得到了松弛主问题(Relaxed Mater Problem, RMP)：

[RMP]

\[\min_{y\in\{y_1,y_2,y_3\}}\quad{q}+f^Ty \]

\[s.t.\quad\quad\quad\quad\quad\quad \]

• 一开始RMP会提供一些比较自私的\(y\)（也就是使目标函数更小的\(y\)），然后我们通过迭代不断地去加cut来约束它。

2.3.2 Benders分解流程

对于RMP的目标函数\(q+f^Ty\)，它其实是提供一个Lower Bound(LB)：因为把MP中的约束都拿走了，RMP的寻优空间就会变大，就会导致它可能会找到一个比原来MP更小的目标函数值。所以不断地往RMP中加入cut，它的可行域就会缩小，它的目标函数值会不断变大，LB就会不断上升。

在RMP输入的一个候选的\(\overline{y}\)下，D-SP有Feasible Bounded的时候才能计算出最优对偶解\(\alpha^*\)以及一个有限的最优目标函数值，而这个值等于SP的最优目标函数值，同时还可以求出\(\overline{y}\)下SP的最优解\(x^*\)。由于SP没有松弛，实际上它此时得到了满足原问题约束的一个可行解\((\overline{y},x^*)\)，它的目标函数值（\(c^Tx^*+f^T\overline{y}\)）等于原问题中一个可行解的目标函数值。由于原问题是 \(min\)，最优值一定不会比某个可行解值更大，所以此时求得的这个目标函数值是一个Upper Bound(UB)。在每次拿到可行解值时，我们都会把这个值和以前最好的那个 UB 比较，如果当前可行解值更小，就更新 UB；如果当前更大，就不更新。所以，迭代的过程就是从当前的可行解跳到下一个更好的可行解，从而使得UB不断地下降。

由此，Benders分解的流程如下：

Step 1： 先松弛MP，得到 RMP。

Step 2： 求解 RMP，得到一个当前的 \(y\) 值。

Step 3： 把这个 \(y\)代入 D-SP 中进行检验。

• 如果 D-SP 是 Feasible Bounded，说明这个 $y $对应的SP可行且有有限最优值。此时加入 optimal cut，目的是修正 RMP 对SP最优值的低估。
• 如果 D-SP 是 Feasible Unbounded，说明这个 \(y\) 会导致原始SP不可行。此时加入 feasibility cut，目的是把这种不合适的 \(y\) 排除掉。

Step 4： 将新生成的 cut 加回 RMP，重新求解，得到新的 \(y\)。

Step 5： 重复上述过程，直到不能再生成新的 cut，或者上下界收敛（也就是\(q+f^Ty=\alpha^T(b-By)+f^Ty\)），此时停止。

3. 一个简单的数值算例

[原问题]

\[min\quad{x}_1+x_2+y \]

\[s.t.\quad{x_1+2x_2+y}\geq4 \]

\[x_1,x_2\geq0 \]

\[y\in[1,2,3,4,5] \]

---

松弛MP，进入Interation1

[RMP(Iteration1)]

\[min\quad{y}+q \]

\[s.t.\quad{y}\in[1,2,3,4,5] \]

随便提供一个对原问题可行的\(y\)给SP(这里选\(y=5,f^Ty=1×5=5\))，得到：

[SP(Iteration1)]

\[min\quad{x}_1+x_2+5 \]

\[s.t.\quad{x_1+2x_2}\geq4-5 \]

\[x_1,x_2\geq0 \]

根据SP(Iteration1)求D-SP(Iteration1)：

[D-SP(Iteration1)]

\[max\quad{-\alpha}+5 \]

\[s.t.\quad{\alpha}\leq1 \]

\[\quad\quad2{\alpha}\leq1 \]

\[\quad\quad\quad{\alpha}\geq0 \]

• D-SP(Iteration1) 有界可行

• 解得：\(\alpha^*=0\)，由此得到\(UB=-0+5=5\)

• Optimal cut：\((b-By)^T\alpha^p_i=(b-By)^T\alpha^*=(4-1×5)×0\leq{q}\)

把Optimal cut加入RMP中，进入Interation2。

---

[RMP(Iteration2)]

\[min\quad{y}+q \]

\[s.t.\quad{q}\geq0 \]

\[\quad\quad\quad\quad\quad\quad\quad{y}\in[1,2,3,4,5] \]

• 求解RMP(Iteration2)，得到：\(q=0,y=1，LB=0+1=1\)

把\(y=1(f^Ty=1×1=1)\)代入SP(Iteration2)：

[SP(Interation2)]

\[min\quad{x}_1+x_2+1 \]

\[s.t.\quad{x_1+2x_2}\geq4-1 \]

\[x_1,x_2\geq0 \]

根据SP(Iteration2)求D-SP(Iteration2)：

[D-SP(Interation2)]

\[max\quad{3\alpha}+1 \]

\[s.t.\quad{\alpha}\leq1 \]

\[\quad\quad2{\alpha}\leq1 \]

\[\quad\quad\quad{\alpha}\geq0 \]

• D-SP(Iteration1) 有界可行

• 解得：\(\alpha^*=0.5,UB=3×0.5+1=2.5\)

• Optimal cut：\((4-1×y)×0.5=2-0.5y\leq{q}\)

把Optimal cut加入RMP中，进入Interation3。

---

[RMP(Iteration3)]

\[min\quad{y}+q \]

\[s.t.\quad{q}\geq0 \]

\[\,\,\quad\quad\quad\quad\quad{q}\geq2-0.5y \]

\[\quad\quad\quad\quad\quad\quad\quad{y}\in[1,2,3,4,5] \]

• 解得：\(y^*=1,q^*=1.5,LB=2.5\)

这时，我们发现：\(LB=UB\)，也就是我们得到了最优解：

• 根据 RMP(Iteration3)，\(y^*=1\)

• 求SP(Interation2)，得到：\(x^*_1=0,x^*_2=1.5\)

• 最优解的目标函数值为\(2.5\)

这样，我们就把两个锅的菜合在一起了。