优化算法 | Benders Decomposition: 一份让你满意的【入门-编程实战-深入理解】的笔记

Jack小新 · 发表于 2024-12-31 12:53:30

本帖最后由 Jack小新于 2025-1-21 15:13 编辑

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.分解算法

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

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

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

3.Benders Decomposition的算法原理

考虑如下的MIP模型：

其中 $c \in R^{1 \times m}, f \in R^{1 \times n}$ ,是行向量， $x, y$ 是列向量，且 $x, y$ 是决策变量，其维度分别为 $m \times 1$ 和 $n \times 1$ 。 $A, B$ 是约束矩阵，其维度分别为 $A \in R^{r \times m}, B \in R^{r \times n}$ 。 $b$ 为右端常数，其维度为 $b \in R^{r \times 1}$ 。注意这里的符号， $R^{r \times m}$ 其实就是表示是一个 $r \times m$ 的一个实数矩阵。

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

● x 是连续型(continuous)决策变量，较为简单；

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

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

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

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

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

我们观察到，如果给定 y 的值，假定为 y¯ ，那么剩余部分的模型就可以写成

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

如何给出 $y$ 的值 $\bar{y}$ 呢？我们可以将关于 $x$ 的部分全部忽略掉，变成下面的模型：

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

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

为什么要交互呢？

原因是：

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

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

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

● Benders optimality cut

● Benders feasibility cut

这两种Cut来源于 SP 的对偶，由于 SP 是LP，因此我们可以将 SP 对偶，变成 Dual SP ，其形式如下：

我们可以观察到：

● Dual SP 的可行性不依赖于 y¯ 的取值。

● 如果强对偶成立(即 SP 和 Dual SP 都有可行解)，则 z_SP=z_{Dual SP}。

● 根据弱对偶性，有 z_SP⩾z_{Dual SP}，即 Dual SP 为 SP 提供了一个下界。

这里我们不加证明地给出Benders optimality cut和Benders feasibility cut的具体形式。

① Benders optimality cuts 的具体形式为

其中， $P$ 表示迭代过程中找到的子问题的对偶问题 $Dual SP$ 的极点(extreme point ，其实就是最优解)的集合， $α_{p}^{i}$ 就表示 $Dual SP$ 的一个极点。

② Benders feasibility cuts 的具体形式为

其中， R 表示迭代过程中找到的子问题的对偶问题 Dual SP 的极射线(extreme ray，其实就是无界射线)的集合， α_rⁱ就表示 Dual SP 的一个极射线。

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

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

① Benders Decomposition 的完整步骤如下：我们将原问题拆分为一个主问题（Master problem, MP）和一个子问题（Subproblem problem, MP）。 $MP$ 如下

其中， $q$ 是辅助变量，是为了衡量目前得到的 $Benders optimality cuts$ 和 $Benders feasibility cuts$ 为原始MIP的目标函数中关于 $x$ 的部分，即 $c x$ 的下界。在MIP的最优解中 $c x$ 的取值 $c x^{*}$ ，是等于在Benders分解中，得到最优解时候对应的 $q^{*}$ 的。

这里， $Benders optimality cuts$ 的具体形式为 $\begin{aligned} (α_{p}^{i})^{T} (b - B y) ⩽ q, \forall i \in P \end{aligned}$
其中， $P$ 表示迭代过程中找到的子问题的对偶问题 $Dual SP$ 的极点的集合， $α_{p}^{i}$ 就表示 $Dual SP$ 的一个极点。
$Benders feasibility cuts$ 的具体形式为 $\begin{aligned} (α_{r}^{i})^{T} (b - B y) ⩽ 0, \forall i \in R \end{aligned}$
其中， $R$ 表示迭代过程中找到的子问题的对偶问题 $Dual SP$ 的极射线(extreme ray)的集合， $α_{r}^{i}$ 就表示 $Dual SP$ 的一个极射线。

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

操作方法1：根据 $\bar{y}$ 构建子问题的对偶问题 $Dual SP$

$\begin{aligned} Dual SP: Q (y) = max_{α} & α^{T} (b - B \bar{y}) \\ s.t. & A^{T} α ⩽ c \\ α free \end{aligned}$

$\begin{aligned} Dual SP: Q (y) = max_{α} & α^{T} (b - B \bar{y}) \\ s.t. & A^{T} α ⩽ c \\ α free \end{aligned}$

$\begin{aligned} Dual SP: Q (y) = max_{α} & α^{T} (b - B \bar{y}) \\ s.t. & A^{T} α ⩽ c \\ α free \end{aligned}$

$\begin{aligned} Dual SP: Q (y) = max_{α} & α^{T} (b - B \bar{y}) \\ s.t. & A^{T} α ⩽ c \\ α free \end{aligned}$

$\begin{aligned} Dual SP: Q (y) = max_{α} & α^{T} (b - B \bar{y}) \\ s.t. & A^{T} α ⩽ c \\ α free \end{aligned}$

$\begin{aligned} Dual SP: Q (y) = max_{α} & α^{T} (b - B \bar{y}) \\ s.t. & A^{T} α ⩽ c \\ α free \end{aligned}$

操作方法1：求解 $Dual SP$ ，直接得到 $α$ 。
如果子问题的对偶 $Dual SP$ 存在有界的可行解，则构建Benders optimality cut ： $\begin{aligned} (α)^{T} (b - B y) ⩽ q \end{aligned}$
如果子问题的对偶 $Dual SP$ 无界，则构建Benders feasibility cut $\begin{aligned} (α)^{T} (b - B y) ⩽ 0 \end{aligned}$
操作方法2：直接求解 $SP$ 的原始形式，然后用求解器得到 $SP$ 的约束的对偶变量 $α_{p}$ 或者极射线 $α_{r}$ 。
如果子问题本身 $SP$ 存在有界的可行解，则构建Benders optimality cut:
$\begin{aligned} (α)^{T} (b - B y) ⩽ q \end{aligned}$
如果子问题本身 $SP$ 无可行解，则构建Benders feasibility cut:
$\begin{aligned} (α)^{T} (b - B y) ⩽ 0 \end{aligned}$

\begin{aligned} (α)^{T} (b - B y) ⩽ 0 \end{aligned}

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

\begin{aligned} L B = f y + q \end{aligned}

\begin{aligned} L B = f y + q \end{aligned}

\begin{aligned} L B = f y + q \end{aligned}

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

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

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

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

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}, \dots, F_{10}$ ，每个基金的收益率分别为1%, 2%, ... 10%。每种基金最多购买100元。
问：小王应该如何投资，才能最大化自己的收益。

本问题的最优解为：

小王应该分别购买100元的 $F_{5}, F_{6}, F_{7}, \dots, F_{10}$ ，剩下的400元存入储蓄账户。

4.2 MIP模型

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

● y ：储蓄账户中要储蓄的金额；

● x_i ：购买第 i 个基金的金额。

则该问题可以建模为下面的MIP模型：

这里，我们直接可以调用Gurobi求解上述模型，具体代码如下：

from gurobipy import *

model = Model('Benders decomposition')

""" create decision variables """
y = model.addVar(lb=0, ub=1000, vtype=GRB.INTEGER, name='y')
x = {}

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

""" set objective function """
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)

""" add constraints """
lhs = LinExpr()
lhs.addTerms(1, y)
for i in range(10):
    lhs.addTerms(1, x)
model.addConstr(lhs <= 1000, name='budget')

model.optimize()
print('\n\n\n')
print('Obj:', model.ObjVal)
print('Saving account : ', y.x)
for i in range(10):
    if(x.x > 0):
        print('Fund ID {}: amount: {}'.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元的 $F_{5}, F_{6}, F_{7}, \dots, F_{10}$ ，剩下的400元存入储蓄账户，最终的本息合计为 $1063.0$ 。

4.3 代码实现：逐步迭代版本ecomposition分解求解小例子：模型分解

虽然上面的模型个非常简单，但是我们还是用Benders Decomposition分解算法来求解一遍上面的模型。由于 $y$ 是整数，而 $x$ 是小数。所以可以将 $y$ 的部分作为主问题，将 $x$ 的部分作为子问题。

因此我们将主问题和子问题分解为下面的形式

主问题

写成具体形式即为：

子问题

写成具体形式即为：

注意，这里 $B = (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)^{T}$

这里我们将主问题稍作变化，我们用 $z$ 表示主问题的目标函数，即 $\begin{aligned} z = f y + q \end{aligned}$

因此，主问题可以写为

这个形式和上面的形式是等价的。只是做了一个等价替换而已。

此外，为了对比上面说的两种等价处理方法，我们把子问题的对偶， $Dual SP$ 也写出来。

具体形式为

这里，我参考的视频https://www.youtube.com/watch?app=desktop&v=vQzpydNOWDY中有错误，我进行了修正。
- Ray (Jian) Zhang. Learn Optimization Visually: Benders Decomposition. url: https://www.youtube.com/watch?

4.4 代码实现：逐步迭代版本

初始化全局的界限 $L B = 0$ , $U B = + \infty$ 。

**主问题（Step 0)**：

这里 $y$ 可以给任意一个正值。比如给 $\bar{y} = 1500$ 。当然，主问题的代码如下

from gurobipy import * MP = Model('Benders decomposition-MP') """ create decision variables """ y = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name='y') z = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name='z') MP.setObjective(z, GRB.MAXIMIZE) MP.addConstr(1000 - y >= 0, name='benders feasibility cut iter 1') MP.optimize() print('\n\n\n') print('Obj:', MP.ObjVal) print('z = %4.1f' % (z.x)) print('y = %4.1f' % (y.x))

第1步

**对偶子问题（Step 1）**

我们可以用求解器求解 $Dual SP1$ ,

from gurobipy import * y_bar = 1500 Dual_SP = Model('Dual SP') """ create decision variables """ alpha_0 = Dual_SP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name='alpha_0') alpha = {} for i in range(10): alpha = Dual_SP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name='alpha_' + str(i)) """ create objective """ obj = LinExpr() obj.addTerms(1000-y_bar, alpha_0) for i in range(10): obj.addTerms(100, alpha) Dual_SP.setObjective(obj, GRB.MINIMIZE) """ add constraints 1-10 """ for i in range(10): Dual_SP.addConstr(alpha_0 + alpha >= 1 + 0.01 * (i+1)) Dual_SP.optimize() print('\n\n\n') print('Model status:', Dual_SP.status) if(Dual_SP.status == 2): print('Obj:', Dual_SP.ObjVal) for i in range(10): if(alpha.x > 0): print('{} = {}'.format(alpha.varName, alpha.x))

求解结果为

Infeasible or unbounded model Model status: 4

此时，我们需要得到模型的极射线。

Gurobi可以很容易地得到极射线，需要使用变量的UnbdRay属性。

手册中是这样描述的：

因此，我们打开InfUnbdInfo，获得极射线。

Dual_SP.setParam('InfUnbdInfo', 1) Dual_SP.optimize() print('\n\n\n') print('Model status:', Dual_SP.status) if(Dual_SP.status == 2): print('Obj:', Dual_SP.ObjVal) for i in range(10): if(alpha.x > 0): print('{} = {}'.format(alpha.varName, alpha.x)) else: print('===== Infeasible or Unbounded information =====') print('extreme ray: {} = {}'.format(alpha_0.varName, alpha_0.UnbdRay)) for i in range(10): print('extreme ray: {} = {}'.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} extreme ray: (α_{0}, α) = (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) \end{aligned}$

于是，可以构造一个Benders feasibility cut

$\begin{aligned} (α_{r}^{i})^{T} (b - B y) ⩾ 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}$

$\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 ⩾ 0 \end{aligned}$

将Benders feasibility cut加入MP，更新MP。

**主问题（step 1）**：

求解 $MP1$ 。

from gurobipy import * MP = Model('Benders decomposition-MP') """ create decision variables """ y = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name='y') z = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.INTEGER, name='z') MP.setObjective(z, GRB.MAXIMIZE) MP.addConstr(1000 - y >= 0, name='benders feasibility cut iter 1') MP.optimize() print('\n\n\n') print('Obj:', MP.ObjVal) print('z = %4.1f' % (z.x)) print('y = %4.1f' % (y.x))

结果为

Obj: 1e+30 z = 1000000000000000019884624838656.0 y = 1000

更新全局的 $U B$ ： $\begin{aligned} U B = z = + \infty . \end{aligned}$

第2步

$\bar{y} = 1000$ ，更新对偶子问题的目标函数如下：

**对偶子问题（Step 2）**：

$\begin{aligned} Dual SP2 : min_{α} & α_{0} (1000 - 1000) + \sum_{i = 1}^{10} 100 α_{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} extreme point: (α_{0}, α) = (1.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) \end{aligned}$

更新全局的 $L B$ ： $\begin{aligned} L B = max {L B, f y + Q (y)} = max {0, 1.045 \times 1000 + 0} = 1045. \end{aligned}$

因为 $Dual SP2$ 有有界可行解，因此可以构造Benders optimality cut。

$\begin{aligned} f y + (α_{p}^{i})^{T} (b - B y) ⩾ z \end{aligned}$

即

$\begin{aligned} 1.045 y + 1.1 (1000 - 1 \cdot y) ⩾ z \\ \to & 1100 - 0.055 y ⩾ z \end{aligned}$

将Benders optimality cut加入MP，更新MP。

**主问题（Step 2）**：

求解 $MP2$ 。

MP.addConstr(1000 - y >= 0, name='benders feasibility cut iter 1') MP.addConstr(1100 - 0.055 * y >= z, name='benders optimality cut iter 2')

结果为

Obj: 1100.0 z = 1100.0 y = -0.0

更新全局的 $U B$ ： $\begin{aligned} U B = z = 1100 \end{aligned}$

此时，Gap为：

第3步

根据上文，此时 $\bar{y} = 0$ ，更新对偶子问题的目标函数如下：

**对偶子问题（Step 3）**：

求解 $Dual SP3$ ,

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} extreme point: (α_{0}, α) = (0, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.10) \end{aligned}$

$Dual SP3$ 有有界可行解。

更新全局的 $L B$ ：

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1045, 0 \times 1000 + 1055} = 1055. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1045, 0 \times 1000 + 1055} = 1055. \end{aligned}$

可以构造Benders optimality cut:

$\begin{aligned} f y + (α_{p}^{i})^{T} (b - B y) ⩾ z \end{aligned}$

即

$\begin{aligned} 1055 + 1.045 y ⩾ z \end{aligned}$

将Benders optimality cut加入MP，更新MP。

**主问题（Step 3）**：

求解 $MP3$ .

MP.addConstr(1000 - y >= 0, name='benders feasibility cut iter 1') MP.addConstr(1100 - 0.055 * y >= z, name='benders optimality cut iter 2') MP.addConstr(1055 + 1.045 * y >= z, name='benders optimality cut iter 3')

结果为

Obj: 1097.0 z = 1097.0 y = 41.0

更新全局的 $U B$ ： $\begin{aligned} U B = z = 1097 \end{aligned}$

此时，Gap为：

第4步

根据上文，此时 $\bar{y} = 41$ ，更新对偶子问题的目标函数如下：

**对偶子问题（Step 4）**：

求解 $Dual SP4$ ,

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} extreme point: (α_{0}, α) = (1.01, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09) \end{aligned}$

$Dual SP4$ 有有界可行解。

更新全局的 $L B$ ：

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1055, 1.045 \times 41 + 1013.59} = 1056.435. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1055, 1.045 \times 41 + 1013.59} = 1056.435. \end{aligned}$

可以构造Benders optimality cut:

$\begin{aligned} f y + (α_{p}^{i})^{T} (b - B y) ⩾ z \end{aligned}$

即

$\begin{aligned} 1055 + 0.035 y ⩾ z \end{aligned}$

将Benders optimality cut加入MP，更新MP。

**主问题（Step 4）**：

求解 $MP4$ .

MP.addConstr(1000 - y >= 0, name='benders feasibility cut iter 1') MP.addConstr(1100 - 0.055 * y >= z, name='benders optimality cut iter 2') MP.addConstr(1055 + 1.045 * y >= z, name='benders optimality cut iter 3') MP.addConstr(1055 + 0.035 * y >= z, name='benders optimality cut iter 4')

结果为

Obj: 1072.0 z = 1072.0 y = 500.0

更新全局的 $U B$ ：

$\begin{aligned} U B = z = 1072 \end{aligned}$

此时，Gap为：

第5步

根据上文，此时 $\bar{y} = 500$ ，更新对偶子问题的目标函数如下：

**对偶子问题（Step 5）**：

求解 $Dual SP5$ ,

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} extreme point: (α_{0}, α) = (1.05, 0, 0, 0, 0, 0, 0.01, 0.02, 0.03, 0.04, 0.05) \end{aligned}$

$Dual SP5$ 有有界可行解。

更新全局的 $L B$ ：

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1056.435, 1.045 \times 500 + 540} = 1062.5. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1056.435, 1.045 \times 500 + 540} = 1062.5. \end{aligned}$

可以构造Benders optimality cut:

$\begin{aligned} f y + (α_{p}^{i})^{T} (b - B y) ⩾ z \end{aligned}$

即

$\begin{aligned} 1065 - 0.005 y ⩾ z \end{aligned}$

将Benders optimality cut加入MP，更新MP。

**主问题（Step 5）**：

求解 $MP5$ .

MP.addConstr(1000 - y >= 0, name='benders feasibility cut iter 1') MP.addConstr(1100 - 0.055 * y >= z, name='benders optimality cut iter 2') MP.addConstr(1055 + 1.045 * y >= z, name='benders optimality cut iter 3') MP.addConstr(1055 + 0.035 * y >= z, name='benders optimality cut iter 4') MP.addConstr(1065 - 0.005 * y >= z, name='benders optimality cut iter 5')

结果为

Obj: 1063.0 z = 1063.0 y = 250.0

更新全局的 $U B$ ：

$\begin{aligned} U B = z = 1063 \end{aligned}$

此时，Gap为：

第6步

根据上文，此时 $\bar{y} = 250$ ，更新对偶子问题的目标函数如下：

**对偶子问题（Step 6）**：

求解 $Dual SP6$ ,

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} extreme point: (α_{0}, α) = (1.03, 0, 0, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07) \end{aligned}$

$Dual SP6$ 有有界可行解。

更新全局的 $L B$ ：

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1062.5, 1.045 \times 250 + 800.5} \\ = max {1062.5, 1061.75} = 1062.5. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1062.5, 1.045 \times 250 + 800.5} \\ = max {1062.5, 1061.75} = 1062.5. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + 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} f y + (α_{p}^{i})^{T} (b - B y) ⩾ z \end{aligned}$

即

$\begin{aligned} 1058 + 0.015 y ⩾ z \end{aligned}$

将Benders optimality cut加入MP，更新MP。

**主问题（Step 6）**：

$\begin{aligned} MP6 : max_{y, z} & z \\ s.t. & y \in N \\ 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \end{aligned}$

$\begin{aligned} MP6 : max_{y, z} & z \\ s.t. & y \in N \\ 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \end{aligned}$

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

结果为

Obj: 1063.0 z = 1063.0 y = 350.0

更新全局的 $U B$ ：

$\begin{aligned} U B = z = 1063 \end{aligned}$

此时，Gap为：

$\begin{aligned} Gap = \frac{U B - L B}{L B} = \frac{1063 - 1062.5}{1062.5} = 0.04 % \end{aligned}$

第7步

根据上文，此时 $\bar{y} = 350$ ，更新对偶子问题的目标函数如下：

**对偶子问题（Step 7）**：

求解 $Dual SP7$ ,

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} extreme point: (α_{0}, α) = (1.04, 0, 0, 0, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06) \end{aligned}$

$Dual SP7$ 有有界可行解。

更新全局的 $L B$ ：

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1062.5, 1.045 \times 350 + 697} \\ = max {1062.5, 1062.75} = 1062.75. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1062.5, 1.045 \times 350 + 697} \\ = max {1062.5, 1062.75} = 1062.75. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1062.5, 1.045 \times 350 + 697} \\ = max {1062.5, 1062.75} = 1062.75. \end{aligned}$

可以构造Benders optimality cut:

$\begin{aligned} f y + (α_{p}^{i})^{T} (b - B y) ⩾ z \end{aligned}$

即

$\begin{aligned} 1061 + 0.005 y ⩾ z \end{aligned}$

将Benders optimality cut加入MP，更新MP。

**主问题（Step 7）**：

求解 $MP7$ .

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

结果为

Obj: 1063.0 z = 1063.0 y = 400.0

更新全局的 $U B$ ：

$\begin{aligned} U B = z = 1063 \end{aligned}$

此时，Gap为：

第8步

根据上文，此时 $\bar{y} = 400$ ，更新对偶子问题的目标函数如下：

**对偶子问题（Step 8）**：

求解 $Dual SP8$ ,

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} extreme point: (α_{0}, α) = (1.04, 0, 0, 0, 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06) \end{aligned}$

$Dual SP8$ 有有界可行解。

更新全局的 $L B$ ：

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1062.75, 1.045 \times 400 + 645} \\ = max {1062.75, 1063} = 1063. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1062.75, 1.045 \times 400 + 645} \\ = max {1062.75, 1063} = 1063. \end{aligned}$

$\begin{aligned} L B & = max {L B, f y + Q (y)} \\ = max {1062.75, 1.045 \times 400 + 645} \\ = max {1062.75, 1063} = 1063. \end{aligned}$

此时 $U B$ 和 $L B$ 相遇了， $G a p = 0$ 。

我们从 $Dual SP$ 的约束的对变量即可获得 $x$ 的取值，代码如下：

print('\n\n\n ==== Solution of SP primal ===== ') cnt = 1 for con in Dual_SP.getConstrs(): print('x {} = {}'.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

到这里，这个小例子的逐步迭代我们就展示完成了。可以看到，我们是对子问题进行了对偶，采用了 $Dual SP$ 。
仅从算法实现的角度，实际上我们也可以不对偶子问题，直接采用子问题本身，即 $SP$ ，来实现Benders Decomposition。只需要使用 Constr.Pi 获得约束的对偶变量，然后拼接出 Cuts 即可。因此并没有必要对子问题进行对偶。
但是从理论推导的角度，还是需要进行对偶的。

4.5 Benders Decomposition过程的Bounds更新

根据上述解释，我们得到，在Benders Decomposition的过程中， $U B$ 和 $L B$ 按照下面的表达式进行更新： $\begin{aligned} L B & = max {L B, f y + Q (y)} \\ U B & = min {U B, f y + α^{T} (b - B \bar{y})} \end{aligned}$ .

我们将上述例子的上下界更新的过程可视化出来，方便读者查看。

4.6 Benders Decomposition的深入理解

我们前面讲到：

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

我们来可视化一下上述例子中，Benders Decomposition在迭代的过程中，recourse是如何以cutting plane的形式补救回来的。

我们只需要画出主问题的可行域的变化即可。

第1步迭代

第2步迭代

第3步迭代

第4步迭代

第5步迭代

第6步迭代

第7步迭代

从上述过程可以看到，每一步子问题都可以提供一个 $z$ 的更紧的上界，一般形式即为 $\begin{aligned} f y + (α_{p}^{i})^{T} (b - B y) ⩾ z \end{aligned}$
本问题的具体形式即为：
$\begin{aligned} 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \\ 1061 + 0.005 y ⩾ z & (Cut 7: optimality cut) \end{aligned}$

$\begin{aligned} 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \\ 1061 + 0.005 y ⩾ z & (Cut 7: optimality cut) \end{aligned}$

$\begin{aligned} 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \\ 1061 + 0.005 y ⩾ z & (Cut 7: optimality cut) \end{aligned}$

$\begin{aligned} 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \\ 1061 + 0.005 y ⩾ z & (Cut 7: optimality cut) \end{aligned}$

$\begin{aligned} 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \\ 1061 + 0.005 y ⩾ z & (Cut 7: optimality cut) \end{aligned}$

$\begin{aligned} 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \\ 1061 + 0.005 y ⩾ z & (Cut 7: optimality cut) \end{aligned}$

$\begin{aligned} 1000 - y ⩾ 0 & (Cut 1: feasibility cut) \\ 1100 - 0.055 y ⩾ z & (Cut 2: optimality cut) \\ 1055 + 1.045 y ⩾ z & (Cut 3: optimality cut) \\ 1055 + 0.035 y ⩾ z & (Cut 4: optimality cut) \\ 1065 - 0.005 y ⩾ z & (Cut 5: optimality cut) \\ 1058 + 0.015 y ⩾ z & (Cut 6: optimality cut) \\ 1061 + 0.005 y ⩾ z & (Cut 7: optimality cut) \end{aligned}$
这些界限，或者说Cuts，就像是对 $z$ 的秋后算账。每一步迭代，都把 $z$ 切一刀，拍一巴掌。相当于把前面忽略掉的 $x$ 的影响，以cutting plane的形式加回来。也就相当于recourse。直到最后的Bound收敛为止。

实际上，如果穷举了所有可能的 $y$ 的取值对应的 $dual SP$ 的极点和极射线的集合 $P$ 和 $R$ ，构造相应的Cuts，并将其全部加入 $MP$ ，则可以直接得到原问题MIP的等价形式。但是这样是非常低效以及不可处理的。因为 $P$ 和 $R$ 的规模是指数级增长的。

虽然 $P$ 和 $R$ 的规模是指数级增长的，但是最后在最优解中，真正binding的却只是一部分。如果能够找到最有可能binding的那部分极点和极射线，这样就可以大大提高求解效率，并且能够保证全局最优性。

Benders Decomposition就是通过迭代的方式，向子问题提供固定了的 $y$ ，并通过求解子问题的对偶问题（或者子问题本身），从而找到相应的极点和极射线。这些找到的极点和极射线，就是非常有可能binding的极点和极射线。

我们通过求解最终的 $MP7$ ，获得约束的对偶变量，来查看一下，究竟最后是哪些Benders feasibility cuts和Benders optimality cuts是binding的。

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

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

具体代码如下：

from gurobipy import * MP = Model('Benders decomposition-MP') """ create decision variables """ y = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name='y') z = MP.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name='z') MP.setObjective(z, GRB.MAXIMIZE) MP.addConstr(1000 - y >= 0, name='benders feasibility cut iter 1') MP.addConstr(1100 - 0.055 * y >= z, name='benders optimality cut iter 2') MP.addConstr(1055 + 1.045 * y >= z, name='benders optimality cut iter 3') MP.addConstr(1055 + 0.035 * y >= z, name='benders optimality cut iter 4') MP.addConstr(1065 - 0.005 * y >= z, name='benders optimality cut iter 5') MP.addConstr(1058 + 0.015 * y >= z, name='benders optimality cut iter 6') MP.addConstr(1061 + 0.005 * y >= z, name='benders optimality cut iter 7') MP.optimize() print('\n\n\n') print('Obj:', MP.ObjVal) print('z = %4.1f' % (z.x)) print('y = %4.1f' % (y.x)) print('\n\n\n ******* Binding information ******* ') for con in MP.getConstrs(): con_dual = con.Pi if(abs(con_dual) > 0): print("cons name: %30s | binding istatus: Yes | Dual: %4.2f" % (con.ConstrName, con.Pi)) else: print("cons name: %30s | binding istatus: No | Dual: %4.2f" % (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

提供了非常紧凑的bound，在最优解中处于binding状态。

其余的Cuts (也就是下面的Cuts)

都是冗余的，是中间迭代过程中产生的暂时性有帮助的Cuts。是我们奔向最优解的垫脚石。最后真正的MVP属于Cut 5和Cut 7.

我们可以验证一下：

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

发现解为：

Obj: 1063.0 z = 1063.0 y = 400.0

也就是说，如果我们找的足够准，直接两步就可以找到Cut 5和Cut 7，那我们2步就可以收敛。

我们只画出Cut 5和Cut 7，如下

小结

本文首先综述了求解MIP的常用算法的主要思想以及适用场景。然后简要介绍了Benders Decomposition的原理，并附以直观易懂的解释。

此外，我们以Youtube上一个视频教程的例子为例，详细的展示了Benders Decomposition的迭代过程。并且我们用Python调用Gurobi，实现了整个过程，并提供了非常详尽的解释。

最后我们从偏理论的角度，通过可视化的手段，将这个过程进行了完整的回顾，并给出了Benders Decomposition的较为深入的理解。

本文的介绍是Benders Decomposition算法的全套学习资料，非常适合想要从理论、公式推导、代码实现3个方面达到全方位入门Benders Decomposition的小伙伴们学习钻研。

参考文献

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

本文转载自知乎博主：运小筹