详解 Benders 分解与一个算例的 python 代码

听说过 benders 分解几年了，在生产管理、路径规划与选址问题中经常应用，一直没有细看，最近论文里面也见到，还是有必要了解一下它的基本思想与用法的。

1. 基本思想
2. 线性规划模型与对偶模型
3. 线性规划的几何意义
- 3.1 多面体
- 3.2 极点与极射线
- 3.3 极射线、有界性、可行解的一些性质
- 3.4 Minkowski 定理
4 Bender's 分解中的各种 problem
5. 迭代过程中的 subproblem 与 master problem
6. 迭代步骤
7. 具体例子求解
8. python 代码

1. 基本思想

benders 分解比较适合求解一些大规模的混合整数规划问题，它的基本思想是：将原问题的求解变量分为两部分，同时也将原问题分解成一个两阶段问题：主问题与子问题；随便初始化子问题的目标函数值，然后求解主问题得到第一部分变量的值（或者随便初始化第一部分变量的值），再将第一部分变量的值代入子问题求解第二部分变量的值；根据子问题是否可行解以及无界解的情况，更新约束并添加到主问题中，主问题的上下界也会发生变化；反复迭代，直到主问题的上界与下界十分接近为止。

注：现在 benders 有很多变种了，上面是经典 Benders 的基本思路。

例如下面一个混合整数规划问题：

问题 P1：
max⁡4x1+7x2+2y1−3y2+y3s.t.2x1+4x2+4y1−2y2+3y3≤123x1+5x2+2y1+3y2−y3≤10x1≤2,x2≤2y1≥0,y2≥0,y3≥0x1,x2为整数\begin{aligned} &\max\quad &4x_1+7x_2+2y_1-3y_2+y_3\\ &s.t.&\\ &&2x_1+4x_2+4y_1-2y_2+3y_3\leq 12\\ && 3x_1+5x_2+2y_1+3y_2-y_3\leq 10\\ &&x_1\leq 2, x_2\leq 2\\ &&y_1\geq 0,y_2\geq 0, y_3\geq 0\\ && x_1, x_2 ~为整数 \end{aligned} maxs.t.4x1+7x2+2y1−3y2+y32x1+4x2+4y1−2y2+3y3≤123x1+5x2+2y1+3y2−y3≤10x1≤2,x2≤2y1≥0,y2≥0,y3≥0x1,x2 为整数

在用 benders 求解这个问题之前，先了解一些数学概念。

2. 线性规划模型与对偶模型

一般的混合整数规划问题可以抽象为下面的表达：
min⁡cTx+dTys.t.Ax+By≥by≥0x∈X\begin{aligned} &\min\quad &\bf{c^\mathrm{T}x+d^\mathrm{T}y}\\ &s.t.\\ &&\bf{Ax+By\geq b}\\ && \bf{y\geq 0}\\ && \bf{x\in X} \end{aligned} mins.t.cTx+dTyAx+By≥by≥0x∈X

其中，X\bf{X}X 为变量 x\bf{x}x 的定义域，例如 X\bf{X}X 可以为整数集 Z\mathbb{Z}Z。对于一个固定值 x‾\bf\overline{x}x，剩余问题(residual problem)可以表示为：
min⁡dTy+cTx‾s.t.By≥b−Ax‾y≥0\begin{aligned} &\min\quad &&\bf{d^\mathrm{T}y+c^\mathrm{T}\overline{x}}\\ &s.t.\\ &&&\bf{By\geq b-A\overline{x}}\\ &&& \bf{y\geq 0} \end{aligned} mins.t.dTy+cTxBy≥b−Axy≥0

剩余问题的对偶问题为：
max⁡(b−Ax‾)u+cTx‾s.t.BTu≥du≥0\begin{aligned} &\max\quad &&\bf{(b-A\overline{x})u+c^\mathrm{T}\overline{x}}\\ &s.t.\\ &&&\bf{B^\mathrm{T}u\geq d}\\ &&& \bf{u\geq 0} \end{aligned} maxs.t.(b−Ax)u+cTxBTu≥du≥0

3. 线性规划的几何意义

3.1 多面体

线性规划的定义域 {x∣Ax≥b}\bf \{x|Ax\geq b\}{x∣Ax≥b} 为一个多面体，例如下面的图形为多面体：

从图形中可以看出，多面体是一个凸集。

3.2 极点与极射线

极点

多面体的顶点也被称为极点（extreme point），它的定义为：对于多面体 P\bf PP，若在多面体内找不到任意两个不同于 x\bf xx 的点 y,z\bf y, zy,z 以及一个标量 λ∈[0,1]\lambda\in[0, 1]λ∈[0,1]，使得 x=λy+(1−λ)z\bf{x}=\lambda \bf{y}+(1-\lambda)\bf{z}x=λy+(1−λ)z，则 x\bf xx 为极点。（即 x\bf xx 不在多面体内任何两点的连线里）

极点可以通过求解这个线性规划，得出的最优解即为一个极点。因为线性规划的最优解总存在于极点上。或者令 n 个线性无关的约束条件取等号，解方程得出，其中 n 为多面体的维度。

射线

多面体的射线定义为：对于一个向量 r\bf rr，若对于多面体 P\bf PP 的任何一点 x\bf xx 以及任意非负标量 λ\lambdaλ，都有 x+λr∈P\bf x+\lambda r\in Px+λr∈P，则称 r\bf rr 为多面体 P\bf PP 的射线。

例如，下面这个多面体中的箭头，都是该多面体的射线。

极射线

多面体的极射线定义为：对于多面体 P\bf PP 中任何两个线性无关的射线 r1\bf r_1r1，r2\bf r_2r2，若射线 r\bf rr 不能表示为 r=λ1r1+λ2r2\bf{r}= \lambda_1\bf{r_1}+\lambda_2\bf{r_2}r=λ1r1+λ2r2，其中 λ1,λ2\lambda_1, \lambda_2λ1,λ2 为大于零的标量，则 r\bf rr 为多面体的极射线。

几何意义：多面体边界上的射线是极射线。

要求出极射线，要用到它的另一个定义：对于一个多面体 {x∣Ax≥b}\bf\{x|Ax\geq b\}{x∣Ax≥b}，它的多面锥（又称退化锥，因为是从多面体退化的）为 {x∣Ax≥0}\bf\{x|Ax\geq 0\}{x∣Ax≥0}，假设 x\bf xx 有 nnn 个维度，则 n-1 个线性无关的约束条件取等号，同时满足剩下的那个条件，则可以求出一个极射线。

注：极射线是一个向量，因此极射线可以有多个解，可以通过求解一个带辅助变量的线性规划问题得出（辅助变量为一个0-1变量，目标函数为这个辅助变量值最大，约束条件中：随便选取 n-1 个约束取等号放进去，剩下那个约束也放进去，随便让其中一个变量等于 1）。

根据《Introduction to linear optimization》一书 179 页，存在无界解时，单纯形表中，判断无界解的那一列的向量取负号，即 d=−cBB−Ajd=-c_BB^-A_jd=−cBB−Aj，即为一个极射线。

3.3 极射线、有界性、可行解的一些性质

一些定理或推论在理解 benders 添加 cuts 时会用到：

弱对偶性：若原问题无界，则对偶问题无可行解；若对偶问题无界，则原问题无可行解；若原问题无可行解，则对偶问题无界或无可行解。
一个多面体若存在极射线，则该多面体无界。
对于一个线性规划问题：min⁡{cTx∣Ax≥b}\bf\min\{c^\mathrm{T}x| Ax\geq b\}min{cTx∣Ax≥b}，假设可行域至少存在一个极点，若该线性规划问题无界，则等同于存在一个极射线 d\bf dd，使得 cTd≤0\bf c^Td\leq 0cTd≤0。
同理，对于一个最大化的线性规划问题：max⁡{cTx∣Ax≤b}\bf\max\{c^\mathrm{T}x| Ax\leq b\}max{cTx∣Ax≤b}，假设可行域至少存在一个极点，若该线性规划问题无界，则等同于存在一个极射线 d\bf dd，使得 cTd≥0\bf c^\mathrm{T}d\geq 0cTd≥0。
对于一个线性规划问题：min⁡{cTx∣Ax≥b}\bf\min\{c^\mathrm{T}x| Ax\geq b\}min{cTx∣Ax≥b}，它有界时，存在可行解的充要条件是：它对偶问题的多面体 {u∣ATu≤c}\bf\{u| A^\mathrm{T}u\leq c\}{u∣ATu≤c} 中所有的极射线 rt\bf r_trt： rtTb≤0,t=1,2,…,T\mathbf{r^\mathrm{T}_tb\leq 0}, ~~ t=1,2,\dots, TrtTb≤0, t=1,2,…,T。(假设有 TTT 个极射线，若对偶问题是无界的，就要找到它的一个极射线，加上这个 cut 约束。)

(结论似乎来自于 Farkas’ Lemma，见课本 Laurence A. Wolsey - Integer programming (2021)，236页，但其他资料里没有看到 Farkas’ Lemma 的这个结论)

对于一个线性规划问题：min⁡{cTx∣Ax≥b}\bf\min\{c^\mathrm{T}x| Ax\geq b\}min{cTx∣Ax≥b}，若它的对偶问题存在可行解，它对偶问题的多面体 {u∣ATu≤c}\bf\{u| A^\mathrm{T}u\leq c\}{u∣ATu≤c} 中 SSS 个极点为 usT,t=1,2,…,S\mathbf{u_s^T}, ~~ t=1,2,\dots, SusT, t=1,2,…,S，设原问题最优值为 η\etaη，则满足：
η=max⁡s=1,2,…,SusTb\eta=\max_{s=1,2,\dots, S}\mathbf{u_s^\mathrm{T}}b η=s=1,2,…,SmaxusTb

上面的式子等价于下面的线性规划：
max⁡ηusTb≤ηs=1,2,…,S\begin{aligned} \max\quad &\eta\\ &\mathbf{u_s^\mathrm{T}}b\leq \eta\quad s=1,2,\dots, S \end{aligned} maxηusTb≤ηs=1,2,…,S

这是由于原问题的目标函数值大于等于对偶问题的最优值，而对偶问题的最优值在它其中一个极点上达到。

3.4 Minkowski 定理

与上面两个性质相关的，还有一个 Minkowski 定理，又称作 Resolution 或 Fication 定理：

任何一个多面体 {x∣Ax≥b}\bf\{x| Ax\geq b\}{x∣Ax≥b}，假设它有 SSS 个极点 xs\bf x_sxs， TTT 个极射线 rt\bf r_trt，则该多面体也等价于下面的表达式：

{x∣x=∑s=1Sλsxs+∑t=1Tvtrt,∑s=1Sλs=1,λs≥0,vt≥0,∀s,∀t}\{\mathbf{x| x}=\sum_{s=1}^S \lambda_s\mathbf{x}_s+\sum_{t=1}^T v_t\mathbf{r}_t, \quad\sum_{s=1}^S \lambda_s=1,\quad \lambda_s\geq 0, v_t\geq 0, \forall s, \forall t\} {x∣x=s=1∑Sλsxs+t=1∑Tvtrt,s=1∑Sλs=1,λs≥0,vt≥0,∀s,∀t}

4 Bender’s 分解中的各种 problem

现在再回到第二部分的线性规划模型中，
原问题：
min⁡cTx+dTys.t.Ax+By≥by≥0x∈X\begin{aligned} &\min\quad &\bf{c^\mathrm{T}x+d^\mathrm{T}y}\\ &s.t.\\ &&\bf{Ax+By\geq b}\\ && \bf{y\geq 0}\\ && \bf{x\in X} \end{aligned} mins.t.cTx+dTyAx+By≥by≥0x∈X

改写为两阶段问题，
第一阶段问题：
min⁡cTx+ϕ(x)s.t.x∈X\begin{aligned} &\min\quad &\bf{c^\mathrm{T}x+\phi(x)}\\ &s.t.\\ && \bf{x\in X} \end{aligned} mins.t.cTx+ϕ(x)x∈X

第二阶段问题 ϕ(x)\bf\phi(x)ϕ(x):
min⁡ϕ(x)=dTys.t.By≥b−Ax‾y≥0\begin{aligned} &\min\quad &&\bf\phi(x)=\bf{d^\mathrm{T}y}\\ &s.t.\\ &&&\bf{By\geq b-A\overline{x}}\\ &&& \bf{y\geq 0} \end{aligned} mins.t.ϕ(x)=dTyBy≥b−Axy≥0

它的对偶问题：
max⁡(b−Ax‾)Tus.t.BTu≤du≥0\begin{aligned} &\max\quad &&\bf{(b-A\overline{x})^\mathrm{T}u}\\ &s.t.\\ &&&\bf{B^\mathrm{T}u\leq d}\\ &&& \bf{u\geq 0} \end{aligned} maxs.t.(b−Ax)TuBTu≤du≥0

设第二阶段问题 ϕ(x)\bf\phi(x)ϕ(x) 的最优值为 η\etaη，对偶问题的极射线为 rtT\mathbf{r}^T_trtT，极点为 usT\mathbf{u}^T_susT.

根据第三部分的性质，原线性规划问题等价于 Bender’s master problem：
(主要源于 Minkowski 定理，查阅的相关资料没有怎么细讲这个结论怎么来的)

min⁡Z=cTx+ηrtT(b−Ax)≤0t=1,2,…T(b−Ax)Tus≤ηs=1,2,…Sx∈X(BMP)\begin{aligned} \min\quad&Z=\mathbf{c^\mathrm{T}x+\eta}\\ &\mathbf{r_\mathrm{t}^\mathrm{T}(b-Ax)\leq 0}\quad t=1,2,\dots T\\ &\mathbf{(b-Ax)^\mathrm{T}u_\mathrm{s}\leq \eta}\quad s=1,2,\dots S\\ &\bf x\in X\tag{BMP} \end{aligned} minZ=cTx+ηrtT(b−Ax)≤0t=1,2,…T(b−Ax)Tus≤ηs=1,2,…Sx∈X(BMP)

由于对偶问题的极点或极射线可以有很多，全 master problem 的约束条件也非常多。在实际求解中，这两个约束条件是逐步添加进去的。

feasibility cut:
rtT(b−Ax)≤0t=1,2,…T(1)\mathbf{r_t^\mathrm{T}(b-Ax)}\leq 0\quad t=1,2,\dots T \tag{1} rtT(b−Ax)≤0t=1,2,…T(1)
这个约束条件保证 BMP 问题可行。

optimality cut:
(b−Ax)Tus≤ηs=1,2,…S(2)\mathbf{(b-Ax)^\mathrm{T}u_s\leq \eta}\quad s=1,2,\dots S \tag{2} (b−Ax)Tus≤ηs=1,2,…S(2)
这个约束条件意味着第二阶段的最优值一定大于等于对偶问题所有极点的目标函数值。

5. 迭代过程中的 subproblem 与 master problem

第一阶段问题结合第二阶段问题的对偶问题，可以写成一个 subproblem，命名为 DSP：
max⁡u(b−Ax‾)Tu+cTx‾s.t.BTu≤du≥0(DSP)\begin{aligned} &\max_\mathbf{u}\quad &&\bf{(b-A\overline{x})^\mathrm{T}u+c^\mathrm{T}\overline{x}}\\ &s.t.\\ &&&\bf{B^\mathrm{T}u\leq d}\\\tag{DSP} &&& \bf{u\geq 0} \end{aligned} umaxs.t.(b−Ax)Tu+cTxBTu≤du≥0(DSP)

若代入DSP问题的 x‾\bf\overline{x}x 值为原问题的可行解，则DSP问题的目标函数值为原问题的一个上界 ZubZ^{ub}Zub。

而 Benders 分解在迭代时求解的主问题 master peroblem 的松弛问题如下，命名为 BR：

min⁡xZlb=cTx+ηcutsx∈PX(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=\mathbf{c^\mathrm{T}x+\eta}\\ \tag{BR} &\text{cuts}\\ &\bf x\in P_X \end{aligned} xminZlb=cTx+ηcutsx∈PX(BR)

其中的 cuts 是逐步放进去的，PX\bf P_XPX 为 x\bf xx 的松弛定义域，例如整数定义域可以松弛为实数。由于约束条件并不完整，所以它求解的目标函数值比原函数更小，为原问题的下界 ZlbZ^{lb}Zlb。

由于原问题一般为混合整数规划问题，一般 Benders 分解要结合分支定界来求解。

6. 迭代步骤

Benders 分解更具体的原理是：反复求解问题 BR 与问题 DSP，更新原问题的上下界，直到上下界很小为止；其中，求解 DSP 问题后，根据问题 DSP 解的情形，添加不同的 cuts 进入到问题 BR 中。

将原问题转化为标准形式，不然影响到对偶问题解的符号。
若 DSP 问题无界，添加 feasibility cut 到 BR 中。
若 DSP 问题有可行解，添加 optimality cut 到 BR 中。
若 DSP 问题无可行解，也表明原问题无可行解（原问题也可能无界，但无界就没必要求了，分支定界时终止这个分支）

实际在求解过程中，还可能结合分支定界对整数变量不断分支求解。

对于分支定界时，要检查 DSP 中的 ϕ(x)=(b−Ax‾)Tu\phi(\bf x)=\bf{(b-A\overline{x})^{\mathrm{T}}u}ϕ(x)=(b−Ax)Tu 是否等于 BR 中求得的 η\etaη 值，若相等，表明 DSP 问题对于给定的 x‾\overline{x}x 值是最优的：此时要检查 x‾\overline{x}x 是否为整数解，若不是，则开始分支，若是，则存储解并更新问题的上界；若 ϕ(x)<η\phi(\bf x)<\etaϕ(x)<η，则添加 feasibility cut 到 BR 中继续迭代。

7. 具体例子求解

对于本文最初的问题 P1，应用 Bender’s 分解的详细步骤为：

原问题：

max⁡4x1+7x2+2y1−3y2+y3s.t.2x1+4x2+4y1−2y2+3y3≤123x1+5x2+2y1+3y2−y3≤10x1≤2,x2≤2y1≥0,y2≥0,y3≥0x1,x2为整数\begin{aligned} &\max\quad &4x_1+7x_2+2y_1-3y_2+y_3\\ &s.t.&\\ &&2x_1+4x_2+4y_1-2y_2+3y_3\leq 12\\ && 3x_1+5x_2+2y_1+3y_2-y_3\leq 10\\ &&x_1\leq 2, x_2\leq 2\\ &&y_1\geq 0,y_2\geq 0, y_3\geq 0\\ && x_1, x_2 ~为整数 \end{aligned} maxs.t.4x1+7x2+2y1−3y2+y32x1+4x2+4y1−2y2+3y3≤123x1+5x2+2y1+3y2−y3≤10x1≤2,x2≤2y1≥0,y2≥0,y3≥0x1,x2 为整数

为了方便，改成最小化问题：

min⁡−4x1−7x2−2y1+3y2−y3s.t.−2x1−4x2−4y1+2y2−3y3≥−12−3x1−5x2−2y1−3y2+y3≥−10y1≥0,y2≥0,y3≥0x1,x2为整数且−x1≥−2,−x2≥−2\begin{aligned} &\min\quad &-4x_1-7x_2-2y_1+3y_2-y_3\\ &s.t.&\\ &&-2x_1-4x_2-4y_1+2y_2-3y_3\geq -12\\ && -3x_1-5x_2-2y_1-3y_2+y_3\geq -10\\ &&y_1\geq 0,y_2\geq 0, y_3\geq 0\\ && x_1, x_2 ~为整数且 -x_1\geq -2, -x_2\geq -2 \end{aligned} mins.t.−4x1−7x2−2y1+3y2−y3−2x1−4x2−4y1+2y2−3y3≥−12−3x1−5x2−2y1−3y2+y3≥−10y1≥0,y2≥0,y3≥0x1,x2 为整数且−x1≥−2,−x2≥−2

迭代时主要用到 b,A,B,c,d\bf b,A, B, c, db,A,B,c,d，几个矩阵的取值分别为：

b=[−12−10]A=[−2−4−3−5]B=[−42−3−2−31]c=[−4−7]d=[−23−1]b=\left[ \begin{aligned} -12\\ -10 \end{aligned} \right] A=\left[ \begin{aligned} -2 ~~&-4\\ -3~~&-5 \end{aligned} \right] B=\left[ \begin{aligned} &-4 &2 &&-3\\ &-2 &-3 &&1 \end{aligned} \right] c=\left[ \begin{aligned} -4\\ -7 \end{aligned} \right] d=\left[ \begin{aligned} -2\\ 3\\ -1 \end{aligned} \right] b=[−12−10]A=[−2 −3 −4−5]B=[−4−22−3−31]c=[−4−7]d=⎣⎢⎡−23−1⎦⎥⎤

它的第二阶段问题 ϕ(x)\phi(\bf x)ϕ(x)为：
max⁡uϕ(x)=(b−Ax‾)Tus.t.BTu≤du≥0\begin{aligned} &\max_\mathbf{u}\quad &&\bf{\phi(x)=(b-A\overline{x})^\mathrm{T}u}\\ &s.t.\\ &&&\bf{B^\mathrm{T}u\leq d}\\ &&& \bf{u\geq 0} \end{aligned} umaxs.t.ϕ(x)=(b−Ax)TuBTu≤du≥0

注意：我们假设它的最优值为 η\etaη，即 η=ϕ(x∗)\eta=\phi(\mathbf{x^\ast})η=ϕ(x∗).

第一步

初始化原问题目标函数值的上界 ZubZ^{ub}Zub 为正无穷，下界 ZlbZ^{lb}Zlb 为负无穷。

由于 master problem 里需要用到 η\etaη (第二阶段问题 ϕ(x)\bf\phi(x)ϕ(x) 的最优值为 η\etaη)，我们随机猜测一个界值。由于第二阶段问题为最大化，随机猜测它的下界为 -200.

master problem BR:
min⁡xZlb=cTx+ηcutsx∈PX(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=\mathbf{c^Tx+\eta}\\ \tag{BR} &\text{cuts}\\ &\bf x\in P_X \end{aligned} xminZlb=cTx+ηcutsx∈PX(BR)

对于这个数例，将 x\bf xx 的整数约束松弛掉，加入 η\etaη 的下界，则 BR 为：

min⁡xZlb=−4x1−7x2+ηη≥−2000≤x1≤2,0≤x2≤2(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=-4x_1-7x_2+\eta\\ \tag{BR} &\eta\geq -200\\ &0\leq x_1\leq 2, 0\leq x_2\leq 2 \end{aligned} xminZlb=−4x1−7x2+ηη≥−2000≤x1≤2,0≤x2≤2(BR)
此时还没有 cuts 约束。求解上面的线性规划问题，得到：x1=2,x2=2,η=−200,Zlb=−222x_1=2, x_2=2,\eta=-200, Z^{lb}=-222x1=2,x2=2,η=−200,Zlb=−222.

（其实也可以随机初始化一组 x\bf xx 值，代入 DSP 问题中求解得到一个初始的 η\etaη 下界）

第二步
求解 subproblem DSP:
max⁡u(b−Ax‾)Tu+cTx‾s.t.BTu≤du≥0(DSP)\begin{aligned} &\max_\mathbf{u}\quad &&\bf{(b-A\overline{x})^\mathrm{T}u+c^\mathrm{T}\overline{x}}\\ &s.t.\\ &&&\bf{B^\mathrm{T}u\leq d}\\\tag{DSP} &&& \bf{u\geq 0} \end{aligned} umaxs.t.(b−Ax)Tu+cTxBTu≤du≥0(DSP)

对于此数例，此时 x‾=(2,2)T\overline{\mathbf{x}}=(2, 2)^\mathrm{T}x=(2,2)T，DSP 模型具体为：
max⁡u0u1+6u2−22s.t.−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)\begin{aligned} &\max_\mathbf{u}\quad &&0 u_1+6u_2-22\\ &s.t.\\ &&&-4u_1-2u_2\leq -2\\\tag{DSP} &&&2u_1-3u_2\leq 3\\ &&&-3u_1+u_2\leq -1\\ &&& u_1\geq 0, u_2\geq0 \end{aligned} umaxs.t.0u1+6u2−22−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)

该 DSP 问题无界，它的一个极射线为 v=(1,3)Tv=(1, 3)^Tv=(1,3)T，就是其中一个约束条件的方向向量。

添加 feasibility cut:
vtT(b−Ax)≤0\mathbf{v_t^\mathrm{T}(b-Ax)}\leq 0 vtT(b−Ax)≤0

即：
11x1+19x2≤4211x_1+19x_2\leq 42 11x1+19x2≤42

添加到 BR 问题中，得到：

min⁡xZlb=−4x1−7x2+ηη≥−20011x1+19x2≤420≤x1≤2,0≤x2≤2(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=-4x_1-7x_2+\eta\\ \tag{BR} &\eta\geq -200\\ &11x_1+19x_2\leq 42\\ &0\leq x_1\leq 2, 0\leq x_2\leq 2 \end{aligned} xminZlb=−4x1−7x2+ηη≥−20011x1+19x2≤420≤x1≤2,0≤x2≤2(BR)
求解该 BR 问题得到：x=(0.36,2)T,η=−200\mathbf{x}=(0.36, 2)^\mathrm{T}, \eta=-200x=(0.36,2)T,η=−200，更新下界为 Zlb=−214.45Z^{lb}=-214.45Zlb=−214.45.

将 x\mathbf{x}x 值代入到 DSP 问题中，
max⁡u−3.27u1+1.09u2−15.45s.t.−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)\begin{aligned} &\max_\mathbf{u}\quad && -3.27 u_1+1.09u_2-15.45\\ &s.t.\\ &&&-4u_1-2u_2\leq -2\\\tag{DSP} &&&2u_1-3u_2\leq 3\\ &&&-3u_1+u_2\leq -1\\ &&& u_1\geq 0, u_2\geq0 \end{aligned} umaxs.t.−3.27u1+1.09u2−15.45−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)

求解得到 u=(0.4,0.2)T\mathbf{u}=(0.4, 0.2)^{\mathrm{T}}u=(0.4,0.2)T，此时 ϕ(x)=−3.27u1+1.09u2=1.09>η=−200\phi(\mathbf{x})= -3.27 u_1+1.09u_2=1.09>\eta=-200ϕ(x)=−3.27u1+1.09u2=1.09>η=−200, 需要添加 cut 返回 BR。

添加 optimality cut:
(b−Ax)Tus≤η\mathbf{(b-Ax)^\mathrm{T}u_s\leq \eta} (b−Ax)Tus≤η

也可以令第一个和第三个约束条件取等号得到极点，optimality cut 具体为：
75x1+135x2−345≤η\frac{7}{5}x_1+\frac{13}{5}x_2-\frac{34}{5}\leq \eta 57x1+513x2−534≤η

添加到 BR 中：
min⁡xZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4275x1+135x2−345≤η0≤x1≤2,0≤x2≤2(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=-4x_1-7x_2+\eta\\ \tag{BR} &\eta\geq -200\\ &11x_1+19x_2\leq 42\\ & \frac{7}{5}x_1+\frac{13}{5}x_2-\frac{34}{5}\leq \eta\\ &0\leq x_1\leq 2, 0\leq x_2\leq 2\\ \end{aligned} xminZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4257x1+513x2−534≤η0≤x1≤2,0≤x2≤2(BR)

得到解 x=(2,1.0526)T\mathbf{x} = (2, 1.0526)^{\mathrm{T}}x=(2,1.0526)T, η=−1.26\eta=-1.26η=−1.26，更新下界 Zlb=−16.4Z^{lb}=-16.4Zlb=−16.4.

代入到DSP 问题中，
max⁡u−3.79u1+1.26u2−15.37s.t.−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)\begin{aligned} &\max_\mathbf{u}\quad && -3.79 u_1+1.26u_2-15.37\\ &s.t.\\ &&&-4u_1-2u_2\leq -2\\\tag{DSP} &&&2u_1-3u_2\leq 3\\ &&&-3u_1+u_2\leq -1\\ &&& u_1\geq 0, u_2\geq0 \end{aligned} umaxs.t.−3.79u1+1.26u2−15.37−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)

得到 u=(0.4,0.2)T\mathbf{u}=(0.4, 0.2)^{\mathrm{T}}u=(0.4,0.2)T，并且 ϕ(x)=−3.79u1+1.26u2=1.26=η\phi(\mathbf{x})= -3.79 u_1+1.26u_2=1.26=\etaϕ(x)=−3.79u1+1.26u2=1.26=η

第三步

由于 x\bf xx 不为整数值，并且 ϕ(x)=η\phi(\mathbf{x})=\etaϕ(x)=η，对 x2x_2x2 进行分支，

分支 x2≤1x_2\leq 1x2≤1 添加到 BR 中，
min⁡xZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4275x1+135x2−345≤η0≤x1≤2,0≤x2≤2x2≤1(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=-4x_1-7x_2+\eta\\ \tag{BR} &\eta\geq -200\\ &11x_1+19x_2\leq 42\\ & \frac{7}{5}x_1+\frac{13}{5}x_2-\frac{34}{5}\leq \eta\\ &0\leq x_1\leq 2, 0\leq x_2\leq 2\\ &x_2\leq 1 \end{aligned} xminZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4257x1+513x2−534≤η0≤x1≤2,0≤x2≤2x2≤1(BR)
得到 x=(2,1)T\mathbf{x} = (2, 1)^{\mathrm{T}}x=(2,1)T，η=−1.4\eta=-1.4η=−1.4。

代入到 DSP 问题中，
max⁡u−4u1+u2−15s.t.−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)\begin{aligned} &\max_\mathbf{u}\quad && -4 u_1+u_2-15\\ &s.t.\\ &&&-4u_1-2u_2\leq -2\\\tag{DSP} &&&2u_1-3u_2\leq 3\\ &&&-3u_1+u_2\leq -1\\ &&& u_1\geq 0, u_2\geq0 \end{aligned} umaxs.t.−4u1+u2−15−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)

得到 u=(0.4,0.2)T\mathbf{u}=(0.4, 0.2)^{\mathrm{T}}u=(0.4,0.2)T，并且 ϕ(x)=−4u1+u2=−1.4=η\phi(\mathbf{x})= -4 u_1+u_2=-1.4=\etaϕ(x)=−4u1+u2=−1.4=η 。由于x\bf xx 值可行，可以更新问题的上界，为 DSP 问题的目标函数值，即 Zub=−16.4Z^{ub}=-16.4Zub=−16.4。发现它与下界 ZlbZ^{lb}Zlb 一样，这个分支终止，并存储该解。

（感觉由于上下界已经一样，下面的步骤可以不要了）

第四步
分支 x2≥2x_2\geq 2x2≥2 添加到 BR 中，
min⁡xZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4275x1+135x2−345≤η0≤x1≤2,0≤x2≤2x2≥2(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=-4x_1-7x_2+\eta\\ \tag{BR} &\eta\geq -200\\ &11x_1+19x_2\leq 42\\ & \frac{7}{5}x_1+\frac{13}{5}x_2-\frac{34}{5}\leq \eta\\ &0\leq x_1\leq 2, 0\leq x_2\leq 2\\ &x_2\geq 2 \end{aligned} xminZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4257x1+513x2−534≤η0≤x1≤2,0≤x2≤2x2≥2(BR)
得到 x=(0.36,2)T\mathbf{x} = (0.36, 2)^{\mathrm{T}}x=(0.36,2)T，η=−1.09\eta=-1.09η=−1.09， Zlb=−16.54Z^{lb}=-16.54Zlb=−16.54。还没有之前的下界好，这个下界不存储。

代入到 DSP 问题中：
max⁡u−3.28u1+1.08u2−15.44s.t.−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)\begin{aligned} &\max_\mathbf{u}\quad && -3.28 u_1+1.08u_2-15.44\\ &s.t.\\ &&&-4u_1-2u_2\leq -2\\\tag{DSP} &&&2u_1-3u_2\leq 3\\ &&&-3u_1+u_2\leq -1\\ &&& u_1\geq 0, u_2\geq0 \end{aligned} umaxs.t.−3.28u1+1.08u2−15.44−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)

得到极点 u=(0.4,0.2)T\mathbf{u}=(0.4, 0.2)^{\mathrm{T}}u=(0.4,0.2)T，并且 ϕ(x)=−3.27u1+1.09u2=−1.09=η\phi(\mathbf{x})= -3.27 u_1+1.09u_2=-1.09=\etaϕ(x)=−3.27u1+1.09u2=−1.09=η，但是 x\bf xx 不是整数解，需要继续分支。

第五步

分支 x1≤0x_1\leq 0x1≤0 添加到 BR 中，
min⁡xZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4275x1+135x2−345≤η0≤x1≤2,0≤x2≤2x2≥2x1≤0(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=-4x_1-7x_2+\eta\\ \tag{BR} &\eta\geq -200\\ &11x_1+19x_2\leq 42\\ & \frac{7}{5}x_1+\frac{13}{5}x_2-\frac{34}{5}\leq \eta\\ &0\leq x_1\leq 2, 0\leq x_2\leq 2\\ &x_2\geq 2\\ &x_1\leq 0 \end{aligned} xminZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4257x1+513x2−534≤η0≤x1≤2,0≤x2≤2x2≥2x1≤0(BR)
得到 x=(0,2)T\mathbf{x} = (0, 2)^{\mathrm{T}}x=(0,2)T，η=−1.6\eta=-1.6η=−1.6， Zlb=−15.6Z^{lb}=-15.6Zlb=−15.6。

代入到 DSP 问题中：
max⁡u−4u1+0u2−14s.t.−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)\begin{aligned} &\max_\mathbf{u}\quad && -4 u_1+0u_2-14\\ &s.t.\\ &&&-4u_1-2u_2\leq -2\\\tag{DSP} &&&2u_1-3u_2\leq 3\\ &&&-3u_1+u_2\leq -1\\ &&& u_1\geq 0, u_2\geq0 \end{aligned} umaxs.t.−4u1+0u2−14−4u1−2u2≤−22u1−3u2≤3−3u1+u2≤−1u1≥0,u2≥0(DSP)

得到极点 u=(0.4,0.2)T\mathbf{u}=(0.4, 0.2)^{\mathrm{T}}u=(0.4,0.2)T，并且 ϕ(x)=−1.6=η\phi(\mathbf{x})= -1.6=\etaϕ(x)=−1.6=η。存储该解，得到的上界为 -15.6，但是这个上界没有之前的好。

第六步
分支 x1≥1x_1\geq 1x1≥1 添加到 BR 中，
min⁡xZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4275x1+135x2−345≤η0≤x1≤2,0≤x2≤2x2≥2x1≥1(BR)\begin{aligned} \min_\mathbf{x}\quad&Z^{lb}=-4x_1-7x_2+\eta\\ \tag{BR} &\eta\geq -200\\ &11x_1+19x_2\leq 42\\ & \frac{7}{5}x_1+\frac{13}{5}x_2-\frac{34}{5}\leq \eta\\ &0\leq x_1\leq 2, 0\leq x_2\leq 2\\ &x_2\geq 2\\ &x_1\geq 1 \end{aligned} xminZlb=−4x1−7x2+ηη≥−20011x1+19x2≤4257x1+513x2−534≤η0≤x1≤2,0≤x2≤2x2≥2x1≥1(BR)

得不到可行解，终止。

综上，最优解为 x∗=(2,1)T\mathbf{x}^\ast = (2, 1)^{\mathrm{T}}x∗=(2,1)T，还可以求出 y∗=(0.1,0,1.2)T\mathbf{y}^\ast = (0.1, 0, 1.2)^{\mathrm{T}}y∗=(0.1,0,1.2)T，最优值 -16.4.

分支定界过程：

上面图片里的目标函数值改为负号。

8. python 代码

上述算例的 python 代码为：

该代码不完善，分支定界没有完全实现，以后再补充。

# -*- coding: utf-8 -*-
"""
Created on Sat Dec 18 20:47:40 2021@author: zhen chen
@Email: chen.zhen5526@gmail.comMIT Licence.Python version: 3.8Description: use benders method for a numerical example:max  4x_1+7x_2+2y_1-3y_2+y_3
s.t.2x_1+4x_2+4y_1-2y_2+3y_3 <= 123x_1+5x_2+2y_1+3y_2-y_3 <= 10x_1 <= 2, x_2 <= 2 y_1 >= 0, y_2 >= 0, y_3 >= 0x_1, x_2 为整数"""from gurobipy import *
import numpy as np
import math# \dsp(x) of the second stage problem, where x is given
def DSP(x):global b global Aglobal Bglobal cglobal dn = len(b) # the number of decision varibales in the second stagetry:m2 = Model("second-stage")u = m2.addMVar(n, vtype=GRB.CONTINUOUS) # gurobi currently support only for 1-D MVar objects, not 2-D# Set objectivecoe = b - A@xm2.setObjective(coe @ u, GRB.MAXIMIZE)m2.addConstr(B.T @ u <= d)m2.update()#m2.setParam('DualReductions', 0)m2.Params.InfUnbdInfo = 1  # or m2.setParam('InfUnbdInfo', 1), Determines whether simplex (and crossover) will compute additional information when a model# is determined to be infeasible or unboundedm2.optimize()if m2.Status == 5: # model is unboundedv = m2.unbdrayreturn m2.Status, v # 若无界返回求解状态与基射线if m2.Status == 2: # model is optimalnita = m2.objValu_value = [u.X[i] for i in range(n)]return m2.Status, u_value, nita  # 若有可行解，返回求解状态与可行值# if m2.Status == 3:#     print("Model is infeasible")#     m2.computeIIS()#     m2.write("model.ilp")except GurobiError as e:print('Error code ' + str(e.errno) + ": " + str(e))b = np.array([-12, -10])
A = np.array([[-2, -4], [-3, -5]])
B = np.array([[-4, 2, -3], [-2, -3, 1]])
c = np.array([-4, -7])
d = np.array([-2, 3, -1])
DSP([1, 1])try:     z_ub = GRB.INFINITYz_lb = -GRB.INFINITY# initialize a lower bound for second stage problemnita_lb = -200 n = len(c) # the number of decision varibales in the first stagem1 = Model("first-stage")x = m1.addMVar(n, ub = 2, vtype=GRB.CONTINUOUS) # gurobi currently support only for 1-D MVar objects, not 2-Dnita = m1.addMVar(1, lb = -GRB.INFINITY, vtype=GRB.CONTINUOUS)m1.addConstr(nita >= nita_lb)m1.setObjective(c @ x + nita, GRB.MINIMIZE)m1.optimize()x_value = [x.X[i] for i in range(n)]z_lb = m1.objVallast_nita_value = -GRB.INFINITY # nita 的一个初始值，若迭代时 nita 的值不发生变化，则开始分支定界all_integer = [False for i in range(n)] # 判断决策变量是不是整数型while abs(z_ub - z_lb) > 1e-6:result = DSP(x_value)if result[0] == 5: # 无界ray = result[1]m1.addConstr(ray@b - ray@A@x <= 0) # 添加 feasibility cutm1.update()m1.optimize()x_value = [x.X[i] for i in range(n)]z_lb = m1.objVal           print()if result[0] == 2: # 可行解u_value = result[1]nita_value = result[2]if abs(last_nita_value - nita_value) < 1e-6: # 开始分支定界for i in range(n):if abs(x_value[i] - math.floor(x_value[i])) < 1e-6:all_integer[i] = Trueelse:        # 这里面的分支定界目前是不完善的，编程实现分支定界要用到递归，# 即模型不断调用自己m1.addConstr(x[i] <= math.floor(x_value[i]))breakif sum(all_integer) == n:z_ub = c @ x_value + nita_valueif abs(z_ub - z_lb) < 1e-6:print('Obj: %g' % m1.objVal)print('x = ')print(x_value)breakm1.addConstr(b@u_value - u_value@A@x <= nita) # 添加 optimality cutm1.update()            m1.optimize()x_value = [x.X[i] for i in range(n)]z_lb = m1.objVallast_nita_value = nita.X[0]  except GurobiError as e:print('Error code ' + str(e.errno) + ": " + str(e))except AttributeError:print('Encountered an attribute error')

参考文献¹² ³

Laurence A. Wolsey - Integer programming, (2021), Second Edition ↩︎
https://bookdown.org/edxu96/matrixoptim/benders-for-standard-milp.html ↩︎
Bertsimas, Dimitris, and John N. Tsitsiklis. Introduction to linear optimization. Vol. 6. Belmont, MA: Athena Scientific, 1997. ↩︎