优化算法 | Benders Decomposition: 一份让你满意的【入门 ...

rustum · 发表于 2022-5-20 14:18

作者：刘兴禄，清华大学，清华-伯克利深圳学院
初次完稿日期：2020.10.08
编辑：张瑞三，四川大学，硕士在读，邮箱：zrssum3@stu.scu.edu.cn
欢迎关注我们的微信公众号 运小筹
本文讲解非常详细，相信一定会对想要入门Benders Decomposition算法的同学们提供非常有效的帮助。本学习笔记是本人2020年整理的，前些天做了一些补充修改。由于读者群很多读者对此都有需求，所以我觉得是时候拿出来分享给大家了。由于笔记太详细，正文大约有1万多字，不过还是希望大家可以仔细读完，正文并没有很多废话，个人感觉都是比较有用的。嘿嘿;)
<hr/>1.大规模混合整数规划求解算法的综述

大规模混合整数规划(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的拓展算法，请读者自行阅读相关文献。
2.分解算法

上面讲到，分解算法，是基于模型的特点，将大问题拆解成为若干小问题。单独求解这些小问题的难度，远低于直接求解原来的大模型。这些小问题之间需要进行一定的交互，从而保证整分解后的模型能够收敛到原来模型的最优解，从而保证分解方法的全局最优性。
这里需要明确， $分解\ne 拆解$ 。分解的操作，实际上是很有艺术性和理论性的。
分解包括 $拆解+耦合$ 。其中拆解就是将问题拆成若干小问题，耦合就是将各个小问题联系起来，形成交互，从而保证全局最优性。
3.Benders Decomposition的算法原理

考虑如下的MIP模型：
$\\\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}$
其中 $c \in \mathbb{R}^{1\times m}, f \in \mathbb{R}^{1\times n}$ ,是行向量，是列向量，且是决策变量，其维度分别为 $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模型具有下面的特点：

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

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

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

因此我们设想：如果我们将这个捣蛋者  先搞定，剩下的部分中，只有  是决策变量，问题就变成了线性规划(Linear programming, LP)了，这不就简单太多了嘛？
Benders Decomposition就是基于这样的思想。
我们观察到，如果给定  的值，假定为  ，那么剩余部分的模型就可以写成 $\\\begin{aligned} \text{SP:} \quad \underset{x}{\min} \quad &cx \\ \text{s.t.} \quad &Ax = b - B\bar{y} \\ & x\geqslant 0 \end{aligned}$
该模型是LP，求解难度比之前的MIP降低了好几档次。因为MIP是NP-hard问题，而LP是多项式时间可精确求解的，不是NP-hard。该问题一般被称之为子问题（Subproblem problem, SP）。
如何给出  的值  呢？我们可以将关于  的部分全部忽略掉，变成下面的模型： $\\ \begin{aligned} \text{MP}_0: \quad \underset{y}{\min} \quad & fy \\ \text{s.t.} \quad & y \in Y \subseteq \mathbb{R}^n \end{aligned}$
求解  非常容易，求解得到的解，即可作为  传给  。该问题一般被称之为主问题（Masterproblem, MP）。  加了下标0表示初始的  (因为后续迭代过程，MP需要加入cut）。
但是，分别求解模型  和  并不能等同于求解了原MIP。要想达到等价地求解原MIP的目的，还需要  和  之间的交互。
为什么要交互呢？
原因是：

删除了所有关于的约束，相当于抛弃了的影响。但是实际上一定会影响。我们通过一种先试错再补救的过程，先删除所有的的相关部分，构成初始的，通过求解，获得一些信息，这些信息反映了对的影响，我们通过某种方式，将对的影响加回到中，进行补救。通过整个过程的迭代，我们最终就可以达到求解原MIP的目的。
其中，补救的步骤，实际上专业术语叫做recourse，反应在具体形式上，就是以cutting plane的形式，将补救措施加回到中。

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

Benders optimality cut
Benders feasibility cut

这两种Cut来源于的对偶，由于是LP，因此我们可以将对偶，变成，其形式如下： $\\ \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}$
我们可以观察到：

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

这里我们不加证明地给出Benders optimality cut和Benders feasibility cut的具体形式。（为什么不加证明呢？是小编不会嘛？当然不是，是太麻烦了，小编在自己写的教材里面给出了非常详细的论述，大家可以去自行查看。或者可以阅读下面的参考文献。）
① Benders optimality cuts 的具体形式为 $\\ \begin{aligned} (\alpha_p^i)^T(b-By) \leqslant q, \quad \forall i \in \mathcal{P} \end{aligned}\\$

其中，表示迭代过程中找到的子问题的对偶问题的极点(extreme point ，其实就是最优解)的集合，就表示的一个极点。
② Benders feasibility cuts 的具体形式为 $\\\begin{aligned} (\alpha_r^i)^T(b-By) \leqslant 0, \quad \forall i \in \mathcal{R} \end{aligned}$
其中，表示迭代过程中找到的子问题的对偶问题的极射线(extreme ray，其实就是无界射线)的集合，就表示的一个极射线。

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

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

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

这里，  的具体形式为 $\\\begin{aligned} (\alpha_p^i)^T(b-By) \leqslant q, \quad \forall i \in \mathcal{P} \end{aligned}$
其中，  表示迭代过程中找到的子问题的对偶问题  的极点的集合，  就表示  的一个极点。
  的具体形式为 $\\ \begin{aligned} (\alpha_r^i)^T(b-By) \leqslant 0, \quad \forall i \in \mathcal{R} \end{aligned}$
其中，  表示迭代过程中找到的子问题的对偶问题  的极射线(extreme ray)的集合，  就表示  的一个极射线。

② 求解主问题，得到的值。并构建子问题。这里有2种等价的操作。

操作方法1：根据构建子问题的对偶问题 $\\\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}$
优点：得到对偶变量方便，直接就是决策变量的取值。
缺点：需要写出子问题的对偶。

操作方法2：这里，也可以不构建，直接只构建即可。即，只需要构建： $\\\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}$
优点：不用写出子问题的对偶，建模方便。
缺点：CPLEX需要用IloCplex.getDuals()获取约束的对偶变量，需要用IloCplex.dualFarkas()函数获得约束的极射线； Gurobi中需要使用Constr.Pi获取约束的对偶变量，需要用Constr.FarkasDual 获取约束的极射线。其实这个也不算啥缺点。

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

操作方法1：求解  ，直接得到 $\alpha$ 。
如果子问题的对偶  存在有界的可行解，则构建Benders optimality cut ：
如果子问题的对偶  无界，则构建Benders feasibility cut
操作方法2：直接求解  的原始形式，然后用求解器得到  的约束的对偶变量 $\alpha_p$ 或者极射线 $\alpha_r$ 。
如果子问题本身  存在有界的可行解，则构建Benders optimality cut:

如果子问题本身  无可行解，则构建Benders feasibility cut:

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

这里需要说明，全局的  和  的计算方法如下： $\\\begin{aligned} &UB = fy + Q(y) \\ &LB = fy + q \end{aligned}$
这是为什么呢？因为 $fy + Q(y)$ ，是给定了 $y=\bar{y}$ 之后，求解  ，得到  的值 $\bar{x}$ ，因此 $(\bar{x}, \bar{y})$ 是原问题MIP的一个可行解。  问题的任意一个可行解，一定是原问题的  ，因此 $\\ \begin{aligned} &UB = fy + Q(y) \end{aligned}$
。而 $fy + q$ 是忽略了  的影响，只是添加了一部分的
Benders optimality cut 和 Benders feasibility cut，并没有将所有的Cuts都加回来。因此它是原问题的一个松弛版本，因此可以为原问题提供一个下界。因此 $\\ \begin{aligned} &LB = fy + q \end{aligned}$

当时，我们就得到了全局最优解，此时正好满足
$\\\begin{aligned} &fy + q = LB = UB = fy + Q(y) \end{aligned}$
也就是 $\\\begin{aligned} q = Q(y) \end{aligned}$

因此，比较简洁的判断Benders Decomposition是否得到最优解的办法，就是判断主问题中的的  取值，是否和子问题的目标函数  相同。
但是大家需要知道为什么是  。背后的原因，在上面的文字中已经给出了详细的理论解释。
直观上来讲，Benders Decomposition的总体框架如下。

图1：Benders Decomposition的总体框架

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

图2：Benders Decomposition算法伪代码：参考自笔者即将出版的教材《运筹优化常用模型、算法与案例实战——Python+Java实现》第17章

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

4.1 问题描述

文献[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元的，剩下的400元存入储蓄账户。

4.2 MIP模型

我们首先对上述问题进行建模。引入下面两组决策变量：

：储蓄账户中要储蓄的金额；
$x_i$ ：购买第 $i$ 个基金的金额。

则该问题可以建模为下面的MIP模型：
$\\\begin{aligned} \begin{matrix} \underset{x, y}{\max} & 1.045 y & + & 1.01x_1 & + & 1.03x_2 & \cdots & 1.10x_{10} & \\ \text{s.t.} (\text{total budget}) & y &+ & x_1&+ & x_2& \cdots & x_{10}& \leqslant &1000 \\ (\text{budget} F_1) & && x_1& & & & & \leqslant &100 \\ (\text{budget} F_2) & && & & x_2& & & \leqslant &100 \\ & && & & & \cdots && \\ (\text{budget} F_{10}) & && & & & & x_{10}& \leqslant &100 \\ & &&x_1, & &x_2, & & x_{10}& \geqslant &0 \\ & y && & & & & & \geqslant &0 \text{ and integer} \end{matrix} \end{aligned}$ 这里，我们直接可以调用Gurobi求解上述模型，具体代码如下：
# * Liu Xinglu
# * hsinglul@163.com
# * Tsinghua University
# * 2022-04-20

from gurobipy import *

model = Model(&#39;Benders decomposition&#39;)

&#34;&#34;&#34; create decision variables &#34;&#34;&#34;
y = model.addVar(lb=0, ub=1000, vtype=GRB.INTEGER, name=&#39;y&#39;)
x = {}

for i in range(10):
x = model.addVar(lb=0, ub=100, vtype=GRB.CONTINUOUS, name=&#39;x_&#39; + str(i))

&#34;&#34;&#34; set objective function &#34;&#34;&#34;
obj = LinExpr()
obj.addTerms(1.045, y)
for i in range(10):
obj.addTerms(1 + 0.01 * (i + 1), x)
model.setObjective(obj, GRB.MAXIMIZE)

&#34;&#34;&#34; add constraints &#34;&#34;&#34;
lhs = LinExpr()
lhs.addTerms(1, y)
for i in range(10):
lhs.addTerms(1, x)
model.addConstr(lhs <= 1000, name=&#39;budget&#39;)

model.optimize()
print(&#39;\n\n\n&#39;)
print(&#39;Obj:&#39;, model.ObjVal)
print(&#39;Saving account : &#39;, y.x)
for i in range(10):
if(x.x > 0):
      print(&#39;Fund ID {}: amount: {}&#39;.format(i+1, x.x))求解结果如下：
Obj: 1063.0
Saving account :  400.0
Fund ID 5: amount: 100.0
Fund ID 6: amount: 100.0
Fund ID 7: amount: 100.0
Fund ID 8: amount: 100.0
Fund ID 9: amount: 100.0
Fund ID 10: amount: 100.0结果与上面介绍的相同。小王应该分别购买100元的  ，剩下的400元存入储蓄账户，最终的本息合计为 $1063.0$   。
4.3 代码实现：逐步迭代版本ecomposition分解求解小例子：模型分解

虽然上面的模型个非常简单，但是我们还是用Benders Decomposition分解算法来求解一遍上面的模型。由于  是整数，而  是小数。所以可以将  的部分作为主问题，将  的部分作为子问题。
因此我们将主问题和子问题分解为下面的形式
主问题
$\\\begin{aligned} \text{MP}: \quad \underset{y}{\max} \quad & fy + q \\ \text{s.t.} \quad & y \in Y \subseteq \mathbb{R}^n \\ \quad & \text{Benders optimality cuts} \\ \quad & \text{Benders feasibility cuts} \\ \quad & q \text{ free} \end{aligned}$
写成具体形式即为：
$\\ \begin{aligned} \text{MP}: \quad \underset{y, q}{\max} \quad & 1.045y + q \\ \text{s.t.} \quad & y \geqslant 0 \text{ and integer} \\ \quad & \text{Benders optimality cuts} \\ \quad & \text{Benders feasibility cuts} \\ \quad & q \text{ free} \end{aligned}$
子问题
$\\\begin{aligned} \text{SP:} \quad Q(y) = \underset{x}{\max} \quad &cx \\ \text{s.t.} \quad &Ax = b - B\bar{y} \\ & x\geqslant 0 \end{aligned}$
写成具体形式即为： $\\\begin{aligned} \begin{matrix} \text{SP}: \quad \underset{x}{\max} & & 1.01x_1 & + & 1.03x_2 & \cdots & 1.10x_{10} & \\ \text{s.t.} & & x_1&+ & x_2& \cdots & x_{10}& \leqslant &1000 - \bar{y} \\ && x_1& & & & & \leqslant &100 \\ && & & x_2& & & \leqslant &100 \\ & && & & & \cdots && \\ && & & & & x_{10}& \leqslant &100 \\ &&x_1, & &x_2, & & x_{10}& \geqslant &0 \end{matrix} \end{aligned}$
注意，这里 $B = (1, 0, 0, 0, 0, 0, 0,0,0,0,0)^T$
这里我们将主问题稍作变化，我们用  表示主问题的目标函数，即 $\\ \begin{aligned} z = fy + q \end{aligned}$
因此，主问题可以写为 $\\\begin{aligned} \text{MP}: \quad \underset{y}{\max} \quad & z \\ \text{s.t.} \quad & y \geqslant 0 \text{ and integer} \\ \quad & fy + (\alpha_p^i)^T (b-By) \geqslant z , \quad \forall i \in \mathcal{P} \\ \quad & (\alpha_r^i)^T (b-By) \geqslant 0 , \quad \forall i \in \mathcal{R} \end{aligned}$
这个形式和上面的形式是等价的。只是做了一个等价替换而已。
此外，为了对比上面说的两种等价处理方法，我们把子问题的对偶，  也写出来。
$\\\begin{aligned} \text{Dual SP:} \quad Q(y) = \underset{\alpha}{\min} \quad & \alpha^T (b-B\bar{y}) \\ \text{s.t.} \quad &A^T\alpha \leqslant c \\ & \alpha \text{ free} \end{aligned}$
具体形式为
$\\\begin{aligned} \text{Dual SP:}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-\bar{y}) + \sum_{i=1}^{10}100\alpha_i && \\ \text{s.t.} \quad &\alpha_0 + \alpha_i \geqslant 1+0.01\cdot i, && \forall i= 1, \cdots, 10. \\ & \alpha_i \geqslant 0, && \forall i= 1, \cdots, 10. \end{aligned}$

这里，我参考的视频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
4.4 代码实现：逐步迭代版本

初始化全局的界限 $LB=0$ , $UB=+\infty$ 。
**主问题（Step 0）**：
$\\\begin{aligned} \text{MP0}: \quad \underset{y, z}{\max} \quad & z \\ \text{s.t.} \quad & y \geqslant 0 \text{ and integer} \end{aligned}$
这里  可以给任意一个正值。比如给 $\bar{y}=1500$ 。当然，主问题的代码如下
from gurobipy import *

MP = Model(&#39;Benders decomposition-MP&#39;)

&#34;&#34;&#34; create decision variables &#34;&#34;&#34;
y = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name=&#39;y&#39;)
z = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name=&#39;z&#39;)

MP.setObjective(z, GRB.MAXIMIZE)

MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)

MP.optimize()
print(&#39;\n\n\n&#39;)
print(&#39;Obj:&#39;, MP.ObjVal)
print(&#39;z = %4.1f&#39; % (z.x))
print(&#39;y = %4.1f&#39; % (y.x))第1步
**对偶子问题（Step 1）**
$\\ \begin{aligned} \text{Dual SP1}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-1500) + \sum_{i=1}^{10}100\alpha_i && \\ \text{s.t.} \quad &\alpha_0 + \alpha_i \geqslant 1+0.01\cdot i, && \forall i= 1, \cdots, 10. \\ & \alpha_0 \geqslant 1.045, && \\ & \alpha_i \geqslant 0, && \forall i= 1, \cdots, 10. \end{aligned}$
我们可以用求解器求解 $\text{Dual SP1}$ ,
from gurobipy import *

y_bar = 1500

Dual_SP = Model(&#39;Dual SP&#39;)

&#34;&#34;&#34; create decision variables &#34;&#34;&#34;
alpha_0 = Dual_SP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name=&#39;alpha_0&#39;)
alpha = {}
for i in range(10):
alpha = Dual_SP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name=&#39;alpha_&#39; + str(i))

&#34;&#34;&#34; create objective &#34;&#34;&#34;
obj = LinExpr()
obj.addTerms(1000-y_bar, alpha_0)
for i in range(10):
obj.addTerms(100, alpha)

Dual_SP.setObjective(obj, GRB.MINIMIZE)

&#34;&#34;&#34; add constraints 1-10 &#34;&#34;&#34;
for i in range(10):
Dual_SP.addConstr(alpha_0 + alpha >= 1 + 0.01 * (i+1))

Dual_SP.optimize()
print(&#39;\n\n\n&#39;)
print(&#39;Model status:&#39;, Dual_SP.status)
if(Dual_SP.status == 2):
print(&#39;Obj:&#39;, Dual_SP.ObjVal)
for i in range(10):
      if(alpha.x > 0):
         print(&#39;{} = {}&#39;.format(alpha.varName, alpha.x))求解结果为
Infeasible or unbounded model

Model status: 4此时，我们需要得到模型的极射线。
Gurobi可以很容易地得到极射线，需要使用变量的UnbdRay属性。
手册中是这样描述的：

图3：gurobi手册中UnbdRay介绍

因此，我们打开InfUnbdInfo，获得极射线。
Dual_SP.setParam(&#39;InfUnbdInfo&#39;, 1)
Dual_SP.optimize()

print(&#39;\n\n\n&#39;)
print(&#39;Model status:&#39;, Dual_SP.status)
if(Dual_SP.status == 2):
print(&#39;Obj:&#39;, Dual_SP.ObjVal)
for i in range(10):
      if(alpha.x > 0):
         print(&#39;{} = {}&#39;.format(alpha.varName, alpha.x))
else:
print(&#39;===== Infeasible or Unbounded information =====&#39;)
print(&#39;extreme ray: {} = {}&#39;.format(alpha_0.varName, alpha_0.UnbdRay))
for i in range(10):
      print(&#39;extreme ray: {} = {}&#39;.format(alpha.varName, alpha.UnbdRay))  结果如下
Model status: 5
===== Infeasible or Unbounded information =====
extreme ray: pi = 1.0
extreme ray: alpha_0 = 0.0
extreme ray: alpha_1 = 0.0
extreme ray: alpha_2 = 0.0
extreme ray: alpha_3 = 0.0
extreme ray: alpha_4 = 0.0
extreme ray: alpha_5 = 0.0
extreme ray: alpha_6 = 0.0
extreme ray: alpha_7 = 0.0
extreme ray: alpha_8 = 0.0
extreme ray: alpha_9 = 0.0我们得到了极射线： $\\\begin{aligned} \text{extreme ray:} \quad (\alpha_0, \boldsymbol{\alpha})=(1, 0, 0, 0,0,0,0,0,0,0,0) \end{aligned}$
于是，可以构造一个Benders feasibility cut
$\\\begin{aligned} (\alpha_r^i)^T (b-By) \geqslant 0 \end{aligned}$
由于
$\\\begin{aligned} &B = (1, 0, 0, 0, 0, 0, 0,0,0,0,0)^T \\ &b = (1000, 100, 100, 100, 100, 100, 100,100,100,100,100)^T \end{aligned}$
即Benders feasibility cut为
$\\ \begin{aligned} 1000 - y \geqslant 0 \end{aligned}$

将Benders feasibility cut加入MP，更新MP。
**主问题（Step 1）**：

求解 $\text{MP1}$ 。
from gurobipy import *

MP = Model(&#39;Benders decomposition-MP&#39;)

&#34;&#34;&#34; create decision variables &#34;&#34;&#34;
y = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name=&#39;y&#39;)
z = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name=&#39;z&#39;)

MP.setObjective(z, GRB.MAXIMIZE)

MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)

MP.optimize()
print(&#39;\n\n\n&#39;)
print(&#39;Obj:&#39;, MP.ObjVal)
print(&#39;z = %4.1f&#39; % (z.x))
print(&#39;y = %4.1f&#39; % (y.x))结果为
Obj: 1e+30
z = 1000000000000000019884624838656.0
y = 1000更新全局的  ： $\\\begin{aligned} &UB = z = +\infty. \end{aligned}$
第2步

$\bar{y}=1000$ ，更新对偶子问题的目标函数如下：
**对偶子问题（Step 2）**：
$\\\begin{aligned} \text{Dual SP2}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-1000) + \sum_{i=1}^{10}100\alpha_i && \end{aligned}$
求解  ,
y_bar = 1000    结果为
Model status: 2
Obj: 0.0
alpha_0 = 1.1
alpha_0 = 0.0
alpha_1 = 0.0
alpha_2 = 0.0
alpha_3 = 0.0
alpha_4 = 0.0
alpha_5 = 0.0
alpha_6 = 0.0
alpha_7 = 0.0
alpha_8 = 0.0
alpha_9 = 0.0我们得到了极点：
$\\\begin{aligned} \text{extreme point:} \quad (\alpha_0, \boldsymbol{\alpha})=(1.1, 0, 0, 0,0,0,0,0,0,0,0) \end{aligned}$
更新全局的  ： $\\\begin{aligned} &LB = \max\{LB, fy + Q(y)\} = \max\{0, 1.045\times1000+0\}=1045. \end{aligned}$
因为  有有界可行解，因此可以构造Benders optimality cut。

即
$\\\begin{aligned} &\quad 1.045y + 1.1(1000-1\cdot y) \geqslant z \\ \rightarrow & \quad 1100 - 0.055 y \geqslant z \end{aligned}$
将Benders optimality cut加入MP，更新MP。
**主问题（Step 2）**：
$\\ \begin{aligned} \text{MP2}: \quad \underset{y, z}{\max} \quad & z && \\ \text{s.t.} \quad & y \in \mathbb{N} && \\ & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \end{aligned}$
求解 $\text{MP2}$ 。
MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)
MP.addConstr(1100 - 0.055 * y >= z, name=&#39;benders optimality cut iter 2&#39;)结果为
Obj: 1100.0
z = 1100.0
y = -0.0更新全局的  ： $\\\begin{aligned} &UB = z = 1100 \end{aligned}$
此时，Gap为：
$\\\begin{aligned} \text{Gap} = \frac{UB-LB}{LB}=\frac{1100-1045}{1045}=5.26\% \end{aligned}$
第3步
根据上文，此时 $\bar{y}=0$ ，更新对偶子问题的目标函数如下：
**对偶子问题（Step 3）**：
$\\\begin{aligned} \text{Dual SP3}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-0) + \sum_{i=1}^{10}100\alpha_i && \end{aligned}$
求解  ,
y_bar = 0结果为
Model status: 2
Obj: 1055.0
alpha_0 = 0.0
alpha_0 = 1.01
alpha_1 = 1.02
alpha_2 = 1.03
alpha_3 = 1.04
alpha_4 = 1.05
alpha_5 = 1.06
alpha_6 = 1.07
alpha_7 = 1.08
alpha_8 = 1.09
alpha_9 = 1.1我们得到了极点： $\\\begin{aligned} \text{extreme point:} \quad (\alpha_0, \boldsymbol{\alpha})=(0, 1.01, 1.02, 1.03,1.04,1.05,1.06,1.07,1.08,1.09,1.10) \end{aligned}$   有有界可行解。
更新全局的  ： $\\\begin{aligned} LB &= \max\{LB, fy + Q(y)\} \\ &= \max\{1045, 0\times1000+1055\}=1055. \end{aligned}$
可以构造Benders optimality cut:

$\begin{aligned} \quad & fy + (\alpha_p^i)^T (b-By) \geqslant z \end{aligned}$
即
$\\ \begin{aligned} 1055 + 1.045y \geqslant z \end{aligned}$
将Benders optimality cut加入MP，更新MP。
**主问题（Step 3）**：

求解 $\text{MP3}$ .

MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)
MP.addConstr(1100 - 0.055 * y >= z, name=&#39;benders optimality cut iter 2&#39;)
MP.addConstr(1055 + 1.045 * y >= z, name=&#39;benders optimality cut iter 3&#39;)结果为
Obj: 1097.0
z = 1097.0
y = 41.0更新全局的  ： $\\\begin{aligned} &UB = z = 1097 \end{aligned}$
此时，Gap为：
$\\\begin{aligned} \text{Gap} = \frac{UB-LB}{LB}=\frac{1097-1055}{1055}=3.98\% \end{aligned}$
第4步
根据上文，此时 $\bar{y}=41$ ，更新对偶子问题的目标函数如下：
**对偶子问题（Step 4）**：
$\\\begin{aligned} \text{Dual SP4}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-41) + \sum_{i=1}^{10}100\alpha_i && \end{aligned}$
求解  ,
y_bar = 0结果为
Model status: 2
Obj: 1013.59
alpha_0 = 1.01
alpha_0 = 0.0
alpha_1 = 0.010000000000000009
alpha_2 = 0.020000000000000018
alpha_3 = 0.030000000000000027
alpha_4 = 0.040000000000000036
alpha_5 = 0.050000000000000044
alpha_6 = 0.06000000000000005
alpha_7 = 0.07000000000000006
alpha_8 = 0.08000000000000007
alpha_9 = 0.09000000000000008我们得到了极点： $\\\begin{aligned} \text{extreme point:} \quad (\alpha_0, \boldsymbol{\alpha})=(1.01, 0, 0.01, 0.02, 0.03,0.04,0.05,0.06,0.07,0.08,0.09) \end{aligned}$   有有界可行解。
更新全局的  ：
$\\\begin{aligned} LB &= \max\{LB, fy + Q(y)\} \\ &= \max\{1055, 1.045\times 41+1013.59\}=1056.435. \end{aligned}$
可以构造Benders optimality cut:

即
$\\\begin{aligned} 1055 + 0.035y \geqslant z \end{aligned}$
将Benders optimality cut加入MP，更新MP。
**主问题（Step 4）**：
： $\\\begin{aligned} \text{MP4}: \quad \underset{y, z}{\max} \quad & z && \\ \text{s.t.} \quad & y \in \mathbb{N} && \\ & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \\ & 1055 + 1.045y \geqslant z &&\text{ (Cut 3: optimality cut)} \\ & 1055 + 0.035y \geqslant z &&\text{ (Cut 4: optimality cut)} \end{aligned}$
求解 $\text{MP4}$ .
MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)
MP.addConstr(1100 - 0.055 * y >= z, name=&#39;benders optimality cut iter 2&#39;)
MP.addConstr(1055 + 1.045 * y >= z, name=&#39;benders optimality cut iter 3&#39;)
MP.addConstr(1055 + 0.035 * y >= z, name=&#39;benders optimality cut iter 4&#39;)结果为
Obj: 1072.0
z = 1072.0
y = 500.0更新全局的  ：
$\\\begin{aligned} &UB = z = 1072 \end{aligned}$
此时，Gap为：
$\\\begin{aligned} \text{Gap} = \frac{UB-LB}{LB}=\frac{1072-1056.435}{1056.435}=1.47\% \end{aligned}$
第5步
根据上文，此时 $\bar{y}=500$ ，更新对偶子问题的目标函数如下：
**对偶子问题（Step 5）**：
$\\\begin{aligned} \text{Dual SP5}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-500) + \sum_{i=1}^{10}100\alpha_i && \end{aligned}$
求解  ,
y_bar = 500结果为
Model status: 2
Obj: 540.0
alpha_0 = 1.05
alpha_0 = 0.0
alpha_1 = 0.0
alpha_2 = 0.0
alpha_3 = 0.0
alpha_4 = 0.0
alpha_5 = 0.010000000000000009
alpha_6 = 0.020000000000000018
alpha_7 = 0.030000000000000027
alpha_8 = 0.040000000000000036
alpha_9 = 0.050000000000000044我们得到了极点：
$\\ \begin{aligned} \text{extreme point:} \quad (\alpha_0, \boldsymbol{\alpha})=(1.05, 0, 0, 0,0,0,0.01,0.02,0.03,0.04, 0.05) \end{aligned}$

有有界可行解。
更新全局的  ：
$\\ \begin{aligned} LB &= \max\{LB, fy + Q(y)\} \\ &= \max\{1056.435, 1.045\times 500+540\}=1062.5. \end{aligned}$
可以构造Benders optimality cut:

即
$\\ \begin{aligned} 1065 - 0.005y \geqslant z \end{aligned}$
将Benders optimality cut加入MP，更新MP。
**主问题（Step 5）**：
$\\ \begin{aligned} \text{MP5}: \quad \underset{y, z}{\max} \quad & z && \\ \text{s.t.} \quad & y \in \mathbb{N} && \\ & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \\ & 1055 + 1.045y \geqslant z &&\text{ (Cut 3: optimality cut)} \\ & 1055 + 0.035y \geqslant z &&\text{ (Cut 4: optimality cut)} \\ & 1065 - 0.005y \geqslant z &&\text{ (Cut 5: optimality cut)} \end{aligned}$
求解 $\text{MP5}$ .
MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)
MP.addConstr(1100 - 0.055 * y >= z, name=&#39;benders optimality cut iter 2&#39;)
MP.addConstr(1055 + 1.045 * y >= z, name=&#39;benders optimality cut iter 3&#39;)
MP.addConstr(1055 + 0.035 * y >= z, name=&#39;benders optimality cut iter 4&#39;)
MP.addConstr(1065 - 0.005 * y >= z, name=&#39;benders optimality cut iter 5&#39;)结果为
Obj: 1063.0
z = 1063.0
y = 250.0更新全局的  ：
$\\ \begin{aligned} &UB = z = 1063 \end{aligned}$
此时，Gap为：

第6步
根据上文，此时 $\bar{y}=250$ ，更新对偶子问题的目标函数如下：
**对偶子问题（Step 6）**：
$\\\begin{aligned} \text{Dual SP6}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-250) + \sum_{i=1}^{10}100\alpha_i && \end{aligned}$
求解  ,
y_bar = 250结果为
Obj: 800.5
alpha_0 = 1.03
alpha_0 = 0.0
alpha_1 = 0.0
alpha_2 = 0.0
alpha_3 = 0.010000000000000009
alpha_4 = 0.020000000000000018
alpha_5 = 0.030000000000000027
alpha_6 = 0.040000000000000036
alpha_7 = 0.050000000000000044
alpha_8 = 0.06000000000000005
alpha_9 = 0.07000000000000006我们得到了极点： $\\\begin{aligned} \text{extreme point:} \quad (\alpha_0, \boldsymbol{\alpha})=(1.03, 0, 0,0, 0.01,0.02,0.03, 0.04,0.05,0.06,0.07) \end{aligned}$

有有界可行解。
更新全局的  ：
$\\\begin{aligned} LB &= \max\{LB, fy + Q(y)\} \\ &= \max\{1062.5, 1.045\times 250+800.5\} \\ &=\max\{1062.5, 1061.75\}=1062.5. \end{aligned}$
可以构造Benders optimality cut:

即
$\\ \begin{aligned} 1058 + 0.015y \geqslant z \end{aligned}$
将Benders optimality cut加入MP，更新MP。
**主问题（Step 6）**：

求解 $\text{MP6}$ .
MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)
MP.addConstr(1100 - 0.055 * y >= z, name=&#39;benders optimality cut iter 2&#39;)
MP.addConstr(1055 + 1.045 * y >= z, name=&#39;benders optimality cut iter 3&#39;)
MP.addConstr(1055 + 0.035 * y >= z, name=&#39;benders optimality cut iter 4&#39;)
MP.addConstr(1065 - 0.005 * y >= z, name=&#39;benders optimality cut iter 5&#39;)
MP.addConstr(1058 + 0.015 * y >= z, name=&#39;benders optimality cut iter 6&#39;)结果为
Obj: 1063.0
z = 1063.0
y = 350.0更新全局的  ：

此时，Gap为：

第7步
根据上文，此时 $\bar{y}=350$ ，更新对偶子问题的目标函数如下：
**对偶子问题（Step 7）**：
$\\\begin{aligned} \text{Dual SP7}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-350) + \sum_{i=1}^{10}100\alpha_i && \end{aligned}$
求解  ,
y_bar = 350结果为
Model status: 2
Obj: 697.0
alpha_0 = 1.04
alpha_0 = 0.0
alpha_1 = 0.0
alpha_2 = 0.0
alpha_3 = 0.0
alpha_4 = 0.010000000000000009
alpha_5 = 0.020000000000000018
alpha_6 = 0.030000000000000027
alpha_7 = 0.040000000000000036
alpha_8 = 0.050000000000000044
alpha_9 = 0.06000000000000005我们得到了极点：
$\\ \begin{aligned} \text{extreme point:} \quad (\alpha_0, \boldsymbol{\alpha})=(1.04, 0, 0,0, 0, 0.01,0.02,0.03, 0.04,0.05,0.06) \end{aligned}$

有有界可行解。
更新全局的  ：
$\\\begin{aligned} LB &= \max\{LB, fy + Q(y)\} \\ &= \max\{1062.5, 1.045\times 350+697\} \\ &=\max\{1062.5, 1062.75\}=1062.75. \end{aligned}$
可以构造Benders optimality cut:

即
$\\\begin{aligned} 1061 + 0.005y \geqslant z \end{aligned}$
将Benders optimality cut加入MP，更新MP。
**主问题（Step 7）**：
$\\\begin{aligned} \text{MP7}: \quad \underset{y, z}{\max} \quad & z && \\ \text{s.t.} \quad & y \in \mathbb{N} && \\ & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \\ & 1055 + 1.045y \geqslant z &&\text{ (Cut 3: optimality cut)} \\ & 1055 + 0.035y \geqslant z &&\text{ (Cut 4: optimality cut)} \\ & 1065 - 0.005y \geqslant z &&\text{ (Cut 5: optimality cut)} \\ & 1058 + 0.015y \geqslant z &&\text{ (Cut 6: optimality cut)} \\ & 1061 + 0.005y \geqslant z &&\text{ (Cut 7: optimality cut)} \end{aligned}$
求解  .
MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)
MP.addConstr(1100 - 0.055 * y >= z, name=&#39;benders optimality cut iter 2&#39;)
MP.addConstr(1055 + 1.045 * y >= z, name=&#39;benders optimality cut iter 3&#39;)
MP.addConstr(1055 + 0.035 * y >= z, name=&#39;benders optimality cut iter 4&#39;)
MP.addConstr(1065 - 0.005 * y >= z, name=&#39;benders optimality cut iter 5&#39;)
MP.addConstr(1058 + 0.015 * y >= z, name=&#39;benders optimality cut iter 6&#39;)
MP.addConstr(1061 + 0.005 * y >= z, name=&#39;benders optimality cut iter 7&#39;)结果为
Obj: 1063.0
z = 1063.0
y = 400.0更新全局的  ：

此时，Gap为：
$\\\begin{aligned} \text{Gap} = \frac{UB-LB}{LB}=\frac{1063-1062.75}{1062.75}=0.002\% \end{aligned}$
第8步
根据上文，此时 $\bar{y}=400$ ，更新对偶子问题的目标函数如下：
**对偶子问题（Step 8）**：
$\\\begin{aligned} \text{Dual SP8}: \quad \underset{\alpha}{\min} \quad & \alpha_0 (1000-400) + \sum_{i=1}^{10}100\alpha_i && \end{aligned}$
求解  ,
y_bar = 400结果为
Obj: 645.0
alpha_0 = 1.04
alpha_0 = 0.0
alpha_1 = 0.0
alpha_2 = 0.0
alpha_3 = 0.0
alpha_4 = 0.010000000000000009
alpha_5 = 0.020000000000000018
alpha_6 = 0.030000000000000027
alpha_7 = 0.040000000000000036
alpha_8 = 0.050000000000000044
alpha_9 = 0.06000000000000005我们得到了极点： $\\\begin{aligned} \text{extreme point:} \quad (\alpha_0, \boldsymbol{\alpha})=(1.04, 0, 0,0, 0, 0.01,0.02,0.03, 0.04,0.05,0.06) \end{aligned}$

有有界可行解。
更新全局的  ：
$\\\begin{aligned} LB &= \max\{LB, fy + Q(y)\} \\ &= \max\{1062.75, 1.045\times 400+645\} \\ &=\max\{1062.75, 1063\}=1063. \end{aligned}$
此时  和  相遇了， $Gap=0$ 。
我们从  的约束的对变量即可获得  的取值，代码如下：
print(&#39;\n\n\n  ==== Solution of SP primal ===== &#39;)
cnt = 1
for con in Dual_SP.getConstrs():
print(&#39;x {} = {}&#39;.format(cnt, con.Pi))
cnt += 1

  ==== Solution of SP primal =====
x 1 = 0.0
x 2 = 0.0
x 3 = 0.0
x 4 = 0.0
x 5 = 100.0
x 6 = 100.0
x 7 = 100.0
x 8 = 100.0
x 9 = 100.0
x 10 = 100.0因此得到了最优解。最优解为
Obj = 1063
y = 400
x_5 = 100
x_6 = 100
x_7 = 100
x_8 = 100
x_9 = 100
x_10 = 100
到这里，这个小例子的逐步迭代我们就展示完成了。可以看到，我们是对子问题进行了对偶，采用了  。
仅从算法实现的角度，实际上我们也可以不对偶子问题，直接采用子问题本身，即  ，来实现Benders Decomposition。只需要使用 Constr.Pi 获得约束的对偶变量，然后拼接出 Cuts 即可。因此并没有必要对子问题进行对偶。
但是从理论推导的角度，还是需要进行对偶的。
4.5 Benders Decomposition过程的Bounds更新

根据上述解释，我们得到，在Benders Decomposition的过程中，  和  按照下面的表达式进行更新： $\\\begin{aligned} LB &= \max\{LB, fy + Q(y)\} \\ UB &= \min\{UB, fy + \alpha^T(b-B\bar{y}) \} \end{aligned}$ .
我们将上述例子的上下界更新的过程可视化出来，方便读者查看。

图4： Benders分解迭代过程中的上下界的变化

4.6 Benders Decomposition的深入理解

我们前面讲到：

删除了所有关于  的约束，相当于抛弃了  的影响。但是实际上  一定会影响  。我们通过一种先试错，再补救的过程，先删除所有的  的相关部分，构成初始的  ，通过求解  ，获得一些信息，这些信息反映了  对  的影响，我们通过某种方式，将  对  的影响加回到  中，进行补救。通过整个过程的迭代，我们最终就可以达到求解原MIP的目的。其中，补救的步骤，实际上专业术语叫做recourse，反应在具体形式上，就是以cutting plane的形式，将补救措施加回到  中。

我们来可视化一下上述例子中，Benders Decomposition在迭代的过程中，recourse是如何以cutting plane的形式补救回来的。
我们只需要画出主问题的可行域的变化即可。
第1步迭代

图5：第1步迭代

第2步迭代
$\\\begin{aligned} \text{MP2}: \quad \underset{y, z}{\max} \quad & z && \\ \text{s.t.} \quad & y \in \mathbb{N} && \\ & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \end{aligned}$

图6：第2步迭代

第3步迭代

图7：第3步迭代

第4步迭代
$\\ \begin{aligned} \text{MP4}: \quad \underset{y, z}{\max} \quad & z && \\ \text{s.t.} \quad & y \in \mathbb{N} && \\ & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \\ & 1055 + 1.045y \geqslant z &&\text{ (Cut 3: optimality cut)} \\ & 1055 + 0.035y \geqslant z &&\text{ (Cut 4: optimality cut)} \end{aligned}$

图8：第4步迭代

第5步迭代
$\\\begin{aligned} \text{MP5}: \quad \underset{y, z}{\max} \quad & z && \\ \text{s.t.} \quad & y \in \mathbb{N} && \\ & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \\ & 1055 + 1.045y \geqslant z &&\text{ (Cut 3: optimality cut)} \\ & 1055 + 0.035y \geqslant z &&\text{ (Cut 4: optimality cut)} \\ & 1065 - 0.005y \geqslant z &&\text{ (Cut 5: optimality cut)} \end{aligned}$

图9：第5步迭代

第6步迭代

图10：第6步迭代

第7步迭代
$\\ \begin{aligned} \text{MP7}: \quad \underset{y, z}{\max} \quad & z && \\ \text{s.t.} \quad & y \in \mathbb{N} && \\ & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \\ & 1055 + 1.045y \geqslant z &&\text{ (Cut 3: optimality cut)} \\ & 1055 + 0.035y \geqslant z &&\text{ (Cut 4: optimality cut)} \\ & 1065 - 0.005y \geqslant z &&\text{ (Cut 5: optimality cut)} \\ & 1058 + 0.015y \geqslant z &&\text{ (Cut 6: optimality cut)} \\ & 1061 + 0.005y \geqslant z &&\text{ (Cut 7: optimality cut)} \end{aligned}$

图11：第7步迭代

从上述过程可以看到，每一步子问题都可以提供一个  的更紧的上界，一般形式即为
本问题的具体形式即为：
$\\\begin{aligned} & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \\ & 1055 + 1.045y \geqslant z &&\text{ (Cut 3: optimality cut)} \\ & 1055 + 0.035y \geqslant z &&\text{ (Cut 4: optimality cut)} \\ & 1065 - 0.005y \geqslant z &&\text{ (Cut 5: optimality cut)} \\ & 1058 + 0.015y \geqslant z &&\text{ (Cut 6: optimality cut)} \\ & 1061 + 0.005y \geqslant z &&\text{ (Cut 7: optimality cut)} \end{aligned}$
这些界限，或者说Cuts，就像是对  的秋后算账。每一步迭代，都把  切一刀，拍一巴掌。相当于把前面忽略掉的  的影响，以cutting plane的形式加回来。也就相当于recourse。直到最后的Bound收敛为止。
实际上，如果穷举了所有可能的  的取值对应的 $\text{dual SP}$ 的极点和极射线的集合  和  ，构造相应的Cuts，并将其全部加入  ，则可以直接得到原问题MIP的等价形式。但是这样是非常低效以及不可处理的。因为  和  的规模是指数级增长的。
虽然  和  的规模是指数级增长的，但是最后在最优解中，真正binding的却只是一部分。如果能够找到最有可能binding的那部分极点和极射线，这样就可以大大提高求解效率，并且能够保证全局最优性。
Benders Decomposition就是通过迭代的方式，向子问题提供固定了的  ，并通过求解子问题的对偶问题（或者子问题本身），从而找到相应的极点和极射线。这些找到的极点和极射线，就是非常有可能binding的极点和极射线。
我们通过求解最终的  ，获得约束的对偶变量，来查看一下，究竟最后是哪些Benders feasibility cuts和Benders optimality cuts是binding的。

注意，由于MIP是没有dual的，因此我们需要将MP改成线性规划。也就是将  设置成连续变量求解。本问题比较特殊，松弛了整数规划，解释不变的。我这么做纯属是为了查看binding的情况，不代表真实的研究可以这么操作，望读者周知。

判定约束是否 binding ，只需要获得约束的对偶变量，查看其是否为0。如果对偶不为0，说明其 binding，如果为0，实际上是可以认为是非binding的，但是这个不是一定的。我这里就认为是0就为非 binding。
具体代码如下：
from gurobipy import *

MP = Model(&#39;Benders decomposition-MP&#39;)

&#34;&#34;&#34; create decision variables &#34;&#34;&#34;
y = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name=&#39;y&#39;)
z = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name=&#39;z&#39;)

MP.setObjective(z, GRB.MAXIMIZE)

MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)
MP.addConstr(1100 - 0.055 * y >= z, name=&#39;benders optimality cut iter 2&#39;)
MP.addConstr(1055 + 1.045 * y >= z, name=&#39;benders optimality cut iter 3&#39;)
MP.addConstr(1055 + 0.035 * y >= z, name=&#39;benders optimality cut iter 4&#39;)
MP.addConstr(1065 - 0.005 * y >= z, name=&#39;benders optimality cut iter 5&#39;)
MP.addConstr(1058 + 0.015 * y >= z, name=&#39;benders optimality cut iter 6&#39;)
MP.addConstr(1061 + 0.005 * y >= z, name=&#39;benders optimality cut iter 7&#39;)

MP.optimize()
print(&#39;\n\n\n&#39;)
print(&#39;Obj:&#39;, MP.ObjVal)
print(&#39;z = %4.1f&#39; % (z.x))
print(&#39;y = %4.1f&#39; % (y.x))

print(&#39;\n\n\n  ******* Binding information ******* &#39;)
for con in MP.getConstrs():
con_dual = con.Pi
if(abs(con_dual) > 0):
      print(&#34;cons name: %30s  | binding istatus: Yes  | Dual: %4.2f&#34; % (con.ConstrName, con.Pi))
else:
      print(&#34;cons name: %30s  | binding istatus: No | Dual: %4.2f&#34; % (con.ConstrName, con.Pi))结果如下：
Obj: 1063.0
z = 1063.0
y = 400.0

******* Binding information *******
cons name: benders feasibility cut iter 1  | binding istatus: No | Dual: 0.00
cons name:  benders optimality cut iter 2  | binding istatus: No | Dual: 0.00
cons name:  benders optimality cut iter 3  | binding istatus: No | Dual: 0.00
cons name:  benders optimality cut iter 4  | binding istatus: No | Dual: 0.00
cons name:  benders optimality cut iter 5  | binding istatus: Yes  | Dual: -0.50
cons name:  benders optimality cut iter 6  | binding istatus: No | Dual: 0.00
cons name:  benders optimality cut iter 7  | binding istatus: Yes  | Dual: -0.50根据上面的信息，可得：Cut 5和Cut 7是binding的，其余的都没有binding.
也就是Cuts

$\\\begin{aligned} & 1065 - 0.005y \geqslant z &&\text{ (Cut 5: optimality cut)} \\ & 1061 + 0.005y \geqslant z &&\text{ (Cut 7: optimality cut)} \end{aligned}$
提供了非常紧凑的bound，在最优解中处于binding状态。
其余的Cuts (也就是下面的Cuts)

$\\\begin{aligned} & 1000 - y \geqslant 0 &&\text{ (Cut 1: feasibility cut)} \\ & 1100 - 0.055 y \geqslant z &&\text{ (Cut 2: optimality cut)} \\ & 1055 + 1.045y \geqslant z &&\text{ (Cut 3: optimality cut)} \\ & 1055 + 0.035y \geqslant z &&\text{ (Cut 4: optimality cut)} \\ & 1058 + 0.015y \geqslant z &&\text{ (Cut 6: optimality cut)} \end{aligned}$
都是冗余的，是中间迭代过程中产生的暂时性有帮助的Cuts。是我们奔向最优解的垫脚石。最后真正的MVP属于Cut 5和Cut 7.
我们可以验证一下：
# MP.addConstr(1000 - y >= 0, name=&#39;benders feasibility cut iter 1&#39;)
# MP.addConstr(1100 - 0.055 * y >= z, name=&#39;benders optimality cut iter 2&#39;)
# MP.addConstr(1055 + 1.045 * y >= z, name=&#39;benders optimality cut iter 3&#39;)
# MP.addConstr(1055 + 0.035 * y >= z, name=&#39;benders optimality cut iter 4&#39;)
MP.addConstr(1065 - 0.005 * y >= z, name=&#39;benders optimality cut iter 5&#39;)
# MP.addConstr(1058 + 0.015 * y >= z, name=&#39;benders optimality cut iter 6&#39;)
MP.addConstr(1061 + 0.005 * y >= z, name=&#39;benders optimality cut iter 7&#39;)发现解为：
Obj: 1063.0
z = 1063.0
y = 400.0也就是说，如果我们找的足够准，直接两步就可以找到Cut 5和Cut 7，那我们2步就可以收敛。
我们只画出Cut 5和Cut 7，如下

图12：Cut 5和Cut 7

小结

本文首先综述了求解MIP的常用算法的主要思想以及适用场景。然后简要介绍了Benders Decomposition的原理，并附以直观易懂的解释。
此外，我们以Youtube上一个视频教程的例子为例，详细的展示了Benders Decomposition的迭代过程。并且我们用Python调用Gurobi，实现了整个过程，并提供了非常详尽的解释。
最后我们从偏理论的角度，通过可视化的手段，将这个过程进行了完整的回顾，并给出了Benders Decomposition的较为深入的理解。
本文的介绍是Benders Decomposition算法的全套学习资料，非常适合想要从理论、公式推导、代码实现3个方面达到全方位入门Benders Decomposition的小伙伴们学习钻研。
完整集成版本代码获取方式

小编也将本文的例子代码写成了集成的循环构造版本。
需要完整代码的小伙伴，可以通过下面的方式获取完整代码：

代码获取方式：  - 转发该推文至朋友圈，集齐30赞，并将集赞截图发送至公众号后台，即可免费获得。  - 转发该推文至朋友圈，集齐15赞，并打赏10元，之后将集赞截图和打赏截图发送至公众号后台，即可获得。  - 不想集赞的小伙伴，可以直接打赏20元，并将打赏截图发送至公众号后台，即可获得。
注意，小编将在2天之后，统一将代码发送给各位小伙伴。2天之内不发送，望各位读者周知。
<hr/>参考文献

1. Ray (Jian) Zhang. Learn Optimization Visually: Benders Decomposition. url: https://www.youtube.com/watch?app=desktop&v=vQzpydNOWDY
2. Zeki Caner Takin. Benders Decomposition. American Cancer Society, 2011. ISBN 9780470400531.
3. 运筹优化常用模型、算法与案例实战——Python+Java实现. 刘兴禄,  熊望祺, 臧永森, 段宏达, 曾文佳, 陈伟坚. 清华大学出版社(即将出版)
<hr/>往期推文合集

第87篇：优化 | 史上最【皮】数学题求解的新尝试：一种求解器的视角 (Python调用Gurobi实现)
第86篇：优化 | 寻找新的建模视角——从直观解释对偶问题入手：以Cutting Stock Problem的对偶问题为例
第85篇：非线性优化 | 非线性问题线性化以及求解的详细案例及Python+Gurobi求解
第84篇：ORers&#39; Bling Chat | 【高光聊天记录集锦-02】：运小筹读者群里那些热烈的讨论
第83篇：Machine-Learning–Based Column Selection for Column Generation
第82篇：最新！205页运小筹优化理论学习笔记发布(2021年9月--2022年3月)！
第81篇：鲁棒优化】| 补充证明:为什么最优解一定满足取值等于绝对值(论文笔记：The Price of Robustness)
第80篇：【鲁棒优化】| 论文笔记：The Price of Robustness - 列不确定性模型部分的推导和一些思考
第79篇：ORers&#39; Bling Chat | 【高光聊天记录集锦-01】：运小筹读者群里那些热烈的讨论
第78篇：优化| 手把手教你学会杉数求解器(COPT)的安装、配置与测试
第77篇：【教学视频】优化 | 线性化(2)：连续变量 * 0-1变量的线性化
第76篇：【教学视频】优化 | 线性化：两个0-1变量相乘的线性化
第75篇：强化学习实战 | DQN和Double DQN保姆级教程：以Cart-Pole为例
第74篇：强化学习| 概念梳理：强化学习、马尔科夫决策过程与动态规划
第73篇：强化学习实战 | Q-learning求解最短路(SPP)问题
第72篇：鲁棒优化 | 以Coding入门鲁棒优化：以一个例子引入（二）
第71篇：鲁棒优化|基于ROME编程入门鲁棒优化：以一个例子引入(一)
第70篇：优化|含分式的非线性规划求解: 基于Outer Approximation的Branch-and-cut 算法及其Java实现
第69篇：科研工具 | 手把手教你玩转文献管理神器：Endnote
第68篇：相约深圳 | 2022 INFORMS 服务科学国际会议·征稿通知
第67篇：鲁棒优化｜利用rome求解鲁棒优化简单案例：以CVaR不确定性集为例
第66篇：机器学习 | 交通流特征工程小技巧和思考
第65篇：最新！145页运小筹优化理论学习笔记发布(2021年4月--9月)！
第64篇：优化 | 手把手教你用Python调用SCIP求解最优化模型
第63篇：优化 | 随机规划案例：The value of the stochastic solution
第62篇：工具 | draw.io: 科研流程示意图必备大杀器
第61篇：优化 | 开源求解器SCIP的安装与入门教程（Windows+Python）
第60篇：优化|Gurobi处理非线性目标函数和约束的详细案例
第59篇：让出租车更“懂”您的出行
第58篇：测试算例下载：《运筹优化常用模型、算法及案例实战：代码手册》
第57篇：鲁棒优化|分布式鲁棒优化转化为SOCP案例及代码实现(Python+RSOME)
第56篇：鲁棒优化 | 分布式鲁棒优化简单案例及实战(RSOME+Gurobi)
第55篇：【尾款+发货】|【代码手册】《运筹优化常用模型、算法及案例实战：Python+Java
第54篇：深度强化学习之：PPO训练红白机1942
第53篇：简单装配线平衡问题
第52篇：【视频讲解】CPLEX的OPL语言：可视化的优化模型求解解决方案
第51篇：算法 | 基于英雄联盟寻路背景的A星算法及python实现
第50篇：【转发】清华大学深圳国际研究生院2021年物流工程与管理项目优秀大学生夏令营报名通知
第49篇：优化 | 精确算法之分支定界介绍和实现（附Python代码）
第48篇：【顶刊论文速递】综述：服务科学视角下的金融科技
第47篇：【重新发布】|《运筹优化常用模型、算法及案例实战：Python+Java实现》【代码手册】开始预购啦！！！
第46篇：智慧交通可视化：手把手教你用Python做出行数据可视化-osmnx 包入门教程
第45篇：优化 | Pick and delivery problem的介绍与建模实现（二）
第44篇：推一个知乎学弱猹的公众号
第43篇：元启发式算法 | 遗传算法(GA)解决TSP问题(Python实现)
第42篇：优化|视频详解Python调用Gurobi的应用案例：TSP
第41篇：最新！213页运小筹优化理论系列笔记发布！
第40篇：运小筹第四期月刊发布！
第39篇：开源交通仿真平台SimMobility的安装教程
第38篇：浅谈 | P问题、NP问题、NPC问题、NP-Hard问题
第37篇：一份掏心掏肺的信息素养笔记分享
第36篇：强化学习｜Q-learning (王者荣耀视角)
第35篇：优化｜高级建模方法（Gurobi）：线性化表达小技巧
第34篇：深度强化学习介绍 | Deep Q Network——理论与代码实现
第33篇：优化 | 列生成算法及Java调用cplex实现
第32篇：优化 | Pick and delivery problem的简介与建模实现（一）
第31篇：最新！运小筹月刊-1月份发布！
第30篇：“Learn to Improve”（L2I）：ICLR文章分享 | 运用RL解决VRP问题
第29篇：线性规划求解小程序 | 基于Python 的Cplex 求解器图形化界面
第28篇：运筹学与管理科学TOP期刊揭秘 —TR系列
第27篇：Julia安装与配置Jupyter Notebook
第26篇：期刊征文| IEEE TRANSACTIONS应对COVID-19的特刊征文消息速递
第25篇：两阶段随机规划(Two-stage Stochastic Programming)：一个详细的例子
第24篇：最新！运小筹月刊-12月份发布！
第23篇：Python调用Gurobi：Assignment Problem简单案例
第22篇：初识随机规划：一个小小例子
第21篇：机器学习运用到VRP的若干小知识
第20篇：运筹学与管理科学TOP期刊揭秘 —Service Science
第19篇：手把手教你用Python实现Dijkstra算法（伪代码+Python实现）
第18篇：运筹学与管理科学大揭秘—TOP期刊主编及研究方向一览
第17篇：优化 | 手把手教你用Python实现动态规划Labeling算法求解SPPRC问题
第16篇：代码 | 运小筹GitHub项目上线啦，快来标星收藏吧！！
第15篇：最新！运小筹月刊首次发布！
第14篇：优化| 手把手教你用Python实现Dijkstra算法（附Python代码和可视化）
第13篇：优化|手把手教你用Python调用Gurobi求解最短路问题
第12篇：优化| 手把手教你用Java实现单纯形法
第11篇：优化| 手把手教你用Python调用Gurobi求解VRPTW
第10篇：优化 | 两阶段启发式算法求解带时间窗车辆路径问题（附Java代码）
第9篇：Java调用cplex求解运输问题
第8篇：优化 | 如何优雅地写出大规模线性规划的对偶
第7篇：优化 | TSP中两种不同消除子环路的方法及Callback实现（Python调用Gurobi求解）
第6篇：元启发式算法 | 禁忌搜索(Tabu Search)解决TSP问题(Python代码实现)
第5篇：论文代码复现 | 无人机与卡车联合配送(Python+Gurobi)
第4篇：优化|Shortest Path Problem及其对偶问题的一些探讨(附Python调用Gurobi实现)
第3篇：可视化|Python+OpenStreetMap的交通数据可视化(二)：OpenStreetMap+Python画城市交通路网图
第2篇：优化|最速下降法：详细案例+Python实现
第1篇：Python+networkX+OpenStreetMap实现交通数据可视化(一)：用OpenStreetMap下载地图数据
<hr/>欢迎关注我们的微信公众号 运小筹
运小筹公众号致力于分享运筹优化(LP、MIP、NLP、随机规划、鲁棒优化)、凸优化、强化学习等研究领域的内容以及涉及到的算法的代码实现。编程语言和工具包括Java、Python、Matlab、CPLEX、Gurobi、SCIP 等。欢迎广大同行投稿！欢迎大家关注我们的公众号“运小筹”以及粉丝群！
可以私信编辑 MunSum3(知乎id) 拉入读者群
<hr/>作者：刘兴禄，清华大学，清华-伯克利深圳学院，博士在读邮箱：hsinglul@163.com,
编辑：张瑞三，四川大学，硕士在读，邮箱：zrssum3@stu.scu.edu.cn

		自动登录	找回密码
密码			立即注册

优化算法 | Benders Decomposition: 一份让你满意的【入门 ...

本帖子中包含更多资源