Python实现规划求解 python求解整数规划

0 介绍

前面介绍的割平面法和分支定界法都是求解整数规划的常用方法，但是面对大规模整数规划/混合整数规划，往往直接采用割平面法和分支定界法求解是不现实的，这时候就需要对大规模整数规划/混合整数规划问题先进行分解和松弛，然后再进一步采用割平面法和分支定界法来帮助求解。目前我个人总结整数规划问题的分解/松弛的主流的方法有如下三种：
1 Benders decomposition (主要思想是行生成+割平面方法)
2 Dantzig-Wolfe decomposition (主要思想其实就是列生成)
3 Lagrangian decomposition (主要思想是 Lagrangian relaxation)
我们今天主要介绍的是 Benders 分解方法，前两种方法我们已经在之前的笔记中介绍过了。

1 Benders分解的Reformulation

考虑如下的线性混合整数规划问题：
$\underset{x,y}{\max}c^Tx+d^Ty\;\;\;\;(1.1) \\ s.t.\ \ Ax+By\le b \;\;\;\;(1.2) \\ x\in R_{+}^{m},\ y\in X\subseteq Z_{+}^{n}\;\;\;\;(1.3)$

其中 $y$ 是离散变量， $x$

在优化问题的求解中，我们直观可以想到的一个方法是先固定住一部分变量，而优化另外一部分变量。例如在上述优化问题（1.1-1.3）中我们可以先将$x$ 作为常数，只把 $y$ 作为决策变量去优化，然后把$y$ 作为常数，只把 $x$ 作为决策变量。在这样一个基本思想的指引下，我们可以将优化问题（1.1-1.3）拆分为2个优化问题：
$W=\underset{y}{\max}\left\{ d^Ty+\phi \left( y \right) :\,\,y\in X\subseteq Z_{+}^{n} \right\} \;\;\;\;(1.4) \\ \phi \left( y \right) =\underset{x}{\max}\left\{ c^Tx\ :\ Ax\le b-By \right\} \;\;\;\;(1.5)$

熟悉 bilevel optimization 的童鞋可能已经看出来就是一个 bilevel optimization 的问题。实际上求解 bilevel optimization 也经常用到行生成的方法，而行生成实际也是 Benders 分解的基本思想。所以 Benders 分解和bilevel optimization有着很密切的联系。在 Benders 分解中我们把优化问题（1.4）称为主问题，把优化问题（1.5）称为子问题。

好了，让我们回到主题上来。可以看到上面的（1.4）就是只包含$y$ 的优化问题，（1.5）就是只包含$x$ 的优化问题。当然需要注意的是优化问题（1.5）是以$y$ 作为参数的，给一个$y$ 的值通过求解优化问题（1.5）就可以得到与之对应的$\phi(y)$。由于$y$ 的值给的不合适可能会让优化问题（1.5）没有可行解。我们不妨假设 $\phi (y) < \infty$，根据 Farkas lemma，优化问题（1.5）有可行解的充要条件为：
$v^t\left( b-By \right) \ge 0,\ t=1,...,T \;\;\;\;(1.6)$
其中$v^t$的极方向。其中极方向集合为
$U=\left\{ v\in R_{+}^{p}:\ A^Tv\ge c \right\} \;\;\;\;(1.7)$
如果大家对 Farkas lemma 不太熟悉的话，那么直接从线性规划的对偶角度来出发也容易得到，优化问题（1.5）有可行解的充要条件。写出优化问题（1.5）的对偶问题为
$\underset{u}{\min}\left\{ \left( b-By \right) ^Tu\ :\ A^Tu\ge c,\ u\in R_{+}^{p} \right\} \;\;\;\;(1.8)$

我们观察发现对偶问题（1.8）的可行域实际上就是极方向集合（1.7）。这是因为线性规划原问题有可行解的充要条件为对偶问题也有可行解。这样进一步我们可以将式（1.5）中的原问题采用式（1.8）的对偶问题的形式来做等价替换可得：

$\phi(y) = \underset{u}{\min}\left\{ \left( b-By \right) ^Tu\ :\ A^Tu\ge c,\ u\in R_{+}^{p} \right\} \;\;\;\;(1.9)$

这里要说一下，从原问题（1.5）折腾了一大圈到目前的对偶问题（1.9）。难道我们直接求解原问题不好吗？为何要绕个圈子去求解对偶问题（1.9）呢？主要原因有以下两点：

1 我们需要对偶变量$u$的信息。之前也说过了原问题可行不可行的关键就是在于极方向集合（1.7），而极方向实际上就是对偶问题的约束，这样转化到对偶问题之后方便我们在主问题里边利用对偶变量的信息$u$来构造cut，这一点在这里大家有点初步体会，在后边会看得更清楚。

2 前面说过了$\phi(y)$是带有参数$y$的一个函数/优化问题，在原问题中 $y$出现在约束上，而在对偶问题中 $y$出现在了目标函数上。显然出现在目标函数上会让带有参数$y$

我们设对偶问题中约束多面体的极点集合为 $\left( u_1,...,u_S \right)$ ，那么优化问题（1.9）进一步可以等价为：
$\phi \left( x \right) =\underset{s=1,...,S}{\min}u_{s}^{T}\left( b-By \right) \;\;\;\;\left( 1.10 \right) \\ \phi \left( y \right) =\underset{\eta}{\max}\left\{ \eta \,\,:\,\,u_{s}^{T}\left( b-By \right) \ge \eta ,\,\,s=1,..,S,\,\,\eta \in R \right\} \;\;\;\;\left( 1.11 \right)$

我们将上式带入到（1.11）中可得：
$W=\underset{y}{\max}\left\{ d^Ty+\underset{\eta}{\max}\left\{ \eta \,\,:\,\,u_{s}^{T}\left( b-By \right) \ge \eta ,\,\,s=1,..,S,\,\,\eta \in R \right\} :\,\,y\in X\subseteq Z_{+}^{n} \right\} \;\;\;\;\left( 1.12 \right)$

上式中是max套max我们可以把max整体提出来合并可得：
$W=\underset{y,\eta}{\max}d^Ty+\eta \;\;\;\;\left( 1.13 \right) \\ s.t.\ \ v_{s}^{T}\left( b-By \right) \ge 0 ,\ \ s=1,..,S \;\;\;\;\left( 1.14 \right) \\y\in X\subseteq Z_{+}^{n}\;\;\;\;\left( 1.15 \right)$

如前所述，我们并不能保证$\phi(y)$ 一定存在可行解，因此还要将保证 $\phi(y)$存在可行解的充要条件（1.6）加入到上述优化问题中，由此可得：
$W=\underset{y}{\max}d^Ty+\eta \;\;\;\;\left( 1.16 \right) \\ s.t.\ \ v_{s}^{T}\left( b-By \right) \ge 0 ,\ \ s=1,..,S \;\;\;\;\left( 1.17 \right) \\ u^T_{t}\left( b-By \right) \ge \eta,\ t=1,...,T \;\;\;\;(1.18)\\ y\in X\subseteq Z_{+}^{n}\;\;\;\;\left( 1.19 \right)$

由于上述约束（1.17-1.18）的数量都是指数级别的，我们是不可能一下子就把所有的约束全部添加进来求解。在实际Benders分解的算法中，我们是逐步逐步将约束（1.17-1.18）添加进来的。
$W=\underset{y}{\max}d^Ty+\eta \;\;\;\;\left( 1.20 \right) \\ s.t.\ \ v_{s}^{T}\left( b-By \right) \ge 0 ,\ \ s \in P \subseteq 1,...,S \;\;\;\;\left( 1.21 \right) \\ u^T_{t}\left( b-By \right) \ge \eta,\ t\in Q\subseteq 1,...,T \;\;\;\;(1.22) \\ y\in X\subseteq Z_{+}^{n}\;\;\;\;\left( 1.23 \right)$

其中集合 $P,Q$

容易看到一条约束实际上就对应一行，Benders分解也就是行生成（逐步添加行的过程），而约束（1.22）实际上就是一种割平面。因此行生成和割平面与Benders分解有着密不可分的联系。

2 Benders分解算法流程

2.1 算法伪代码

Step 1:初始化 $y^*$，初始化上下界 $UB:=+\infty, LB:=-\infty$

Step 2: 若 $UB-LB>\varepsilon$则进入Step 3，否则直接输出最优解

Step 3: 求解对偶问题（1.8）得到 最优解$u^*$

Step 4: 若$u^*$是 unbounded 则添加约束$ v_{s}^{T}\left( b-By \right) \ge 0$到主问题中

Step 5: 若$u^*$是 bounded 则添加约束$u_{s}^{T}\left( b-By \right) \ge \eta$到主问题中， 更新 $LB:=\max\text{\{}LB,\ d^Ty^*+\left( b-By^* \right) ^Tu^*\}$

Step 6: 求解主问题（1.20-1.23）得到最优解 $y^*, \eta^*$，然后更新$UB := d^Ty^*+\eta^*$

2.2 以一个案例进一步说明算法流程

考虑如下的优化问题：
$\max 1.01x_1+1.02x_2+...+1.10x_{10}+1.045y \\ s.t. \;\; x_1+x_2+...+x_{10}+y\le 1000 \\ x_1 \le 100\\ x_2 \le 100\\ ...\\ x_{10} \le 100 \\ x_1,x_2,...,x_{10},y \ge 0$
套用下面的形式：
$\underset{x,y}{\max}c^Tx+d^Ty\;\;\;\;(2.1) \\ s.t.\ \ Ax+By\le b \;\;\;\;(2.2) \\ x\in R_{+}^{m},\ y\in X\subseteq Z_{+}^{n}\;\;\;\;(2.3)$
可得参数为：
$\mathbf{x}=\left( x_1,x_2,...,x_{10} \right) ^T,\ \mathbf{y}=\left( y \right) \\ \mathbf{c}=\left( 1.01,1.02,...,1.10 \right) ^T,\ \mathbf{d}=\left( 1.045 \right)$
$\mathbf{A}=\left( \begin{matrix} 1& ...& 1\\ 1& ...& 0\\ \vdots& \ddots& \vdots\\ 0& ...& 1\\ \end{matrix} \right) ,\ \mathbf{B}=\left( \begin{array}{c} 1\\ 0\\ \vdots\\ 0\\ \end{array} \right) ,\ \mathbf{b}=\left( \begin{array}{c} 1000\\ 100\\ \vdots\\ 100\\ \end{array} \right)$

Step 1: 令 $y^*=1500$，初始化上下界 $UB=+\infty, LB=-\infty$

算法第一次循环：

Step 2: $UB-LB=+\infty >\varepsilon$因此进入 Step 3

Step 3: 带入到子问题的对偶问题（1.8）中可得：
$\underset{u}{\min}\left( b-By^* \right) ^Tu=-500u_1+100u_2+...+100u_{11} \\ s.t.\ \ u_1+u_2\ge 1.01,\\ u_1+u_3\ge 1.02,\\ ...\\ u_1+u_{10} \ge 1.1,\\ u_1 \ge 1.045,\\ u_1,...,u_{11} \ge 0$

最优解为 $u_1=+\infty, u_2,...u_{11}=0$，目标函数为 $-\infty$。因为子问题的对偶问题为unbounded，所以子问题的原问题是 不可行的，可以得到其极方向为 $v=\left( 1,0,...,0 \right) ^T$接下来到Step 4。

Step 4: 若 $u^*$ 是 unbounded 则添加约束 $v_{s}^{T}\left( b-By \right) \ge 0$到主问题中，带入可得 $1000-y \ge 0$

$\underset{y,\eta}{\max}1.045y + \eta \\ 1000-y \ge 0\\ y \in Z^{+},\;\eta \in R$

由此可得最优目标函数为 $+\infty$, 最优解为$y^*=1000, \eta^*=+\infty$，进一步更新 $UB=+\infty$

算法第二次循环：
Step 2: $UB-LB=+\infty>\varepsilon$因此进入 Step 3

Step 3: 带入 $y^*=1000$到子问题的对偶问题（1.8）中可得：

$\underset{u}{\min}\left( b-By^* \right) ^Tu=100u_2+...+100u_{11} \\ s.t.\ \ u_1+u_2\ge 1.01,\\ u_1+u_3\ge 1.02,\\ ...\\ u_1+u_{10} \ge 1.1,\\ u_1,...,u_{11} \ge 0$

最优解为 $u_1=1.1, u_2,...u_{11}=0$，目标函数为 0 。因为子问题的对偶问题为 bounded，所以子问题的原问题是 可行的，接下来到 Step 5。

Step 5: 若 $u^*$是 bounded 则添加约束 $u_{s}^{T}\left( b-By \right) \ge \eta$到主问题中，带入可得 $1100-y \ge \eta$

$\underset{y,\eta}{\max}1.045y + \eta \\ s.t. \;\;1000-y \ge 0\\ 1100-y \ge \eta\\ y \in Z^{+}\;\eta \in R$

由此可得最优目标函数为 $+\infty$，最优解为 $y^*=0, \eta^*=1100$，进一步更新 $UB=1100,LB = 1045$

算法第三次循环：

Step 2: $UB-LB=55>\varepsilon$

Step 3: 带入 $y^*=0$到子问题的对偶问题（1.8）中可得：
$\underset{u}{\min}\left( b-By^* \right) ^Tu=1000u_1+100u_2+...+100u_{11} \\ s.t.\ \ u_1+u_2\ge 1.01,\\ u_1+u_3\ge 1.02,\\ ...\\ u_1+u_{10} \ge 1.1,\\ u_1,...,u_{11} \ge 0$

最优解为 $u_1=0, u_2 =1.01,u_3=1.02,...u_{11}=1.1$

Step 5: 若$u^*$ 是 bounded 则添加约束$u_{s}^{T}\left( b-By \right) \ge \eta$到主问题中，带入可得 $1055 \ge \eta$
$\underset{y, \eta}{\max}1.045y + \eta \\ s.t. \;\;1000-y \ge 0\\ 1100-y \ge \eta\\ 1055 \ge \eta\\ y \in Z^{+}, \;\eta \in R$

由此可得最优目标函数为 1097.745 最优解为$y^*=41, \eta^*=1054.605$，进一步更新$UB=1097.745,LB = 1045$

重复上述迭代过程，直到满足收敛条件为止。

3 Benders分解代码实现

代码对应上一节的算法流程，首先定义初始的优化问题。

N, M = 10, 1 # 初始化决策变量维数
c, d = np.array([1+0.01*i for i in range(1, N+1)]), 1.045  # 初始化系数
A, B = np.vstack((np.ones((1, N)), np.eye(N))), np.array([1 if i == 0 else 0 for i in range(N+1)]).reshape(N+1,1)
b = np.array([1000 if i == 0 else 100 for i in range(N+1)]).reshape(N+1,1)

进入代码的主循环部分：

ub, lb = np.inf, -np.inf # 初始化上下界
MAX_ITER_TIMES, eps = 10, 0.1 # 初始化最大迭代次数和误差

subproblem = Subproblem(N, M) # 定义子问题
subproblem.add_constrs(A, c)
masterproblem = Master(N, M, d) #定义主问题
masterproblem.set_objective()
y = 1500 # 设置y的初始值

for i in range(MAX_ITER_TIMES):
    if ub - lb <= eps:
        break
    subproblem.set_objective(B, b, y)
    subproblem.solve()
    subproblem.write()
    rays, solution_status = subproblem.get_status()
    
    if solution_status == GRB.Status.UNBOUNDED or solution_status == GRB.Status.INF_OR_UNBD:
        masterproblem.add_cut1(B, b, u = rays) # 添加不可行的cut对应式（1.17）
    if solution_status == GRB.Status.OPTIMAL:
        masterproblem.add_cut2(B, b, u = rays)
        lb = max(lb, subproblem.get_objval() + d*y) # 添加可行的cut对应式（1.18）
 
    masterproblem.solve()
    masterproblem.write()
    y = masterproblem.get_solution()
    ub = masterproblem.get_objval()

以上仅展示了代码的主干部分，完整代码可在GitHub下载：https://github.com/WenYuZhi/EasyIntegerProgramming

参考文献：

【1】孙小玲，李端，整数规划，科学出版社，2010
【2】Laurence A. Wolsey, Integer Programming, Wily, 2021
【3】Benders Decomposition: An Easy Example