【整数规划算法】分支定界法及其Python代码实现

运筹优化博士，只做原创博文。更多关于运筹学，优化理论，数据科学领域的内容，欢迎关注我的知乎账号：https://www.zhihu.com/people/wen-yu-zhi-37

分支定界法(Branch and Bound)是整数规划领域最基本的一个算法。该算法的利用了两大最基本最优性条件和分而治之来解决整数规划问题。本节就从这个两大思路展开分支定界法的由来，并且给出具体的代码实现，帮助大家理解分支定界法。更多关于运筹学，最优化，数据科学，机器学习，强化学习的精彩内容可以关注我的知乎账号：https://www.zhihu.com/people/wen-yu-zhi-37

1 最优性条件与分而治之两大思路的结合

求解一般的最优化问题的基本策略是先建立该问题的最优性条件，然后根据最优性条件利用迭代算法寻找满足该条件的可行解。这种基本策略已经在线性规划，二次规划，非线性规划问题中得到了很多使用。而在整数规划问题中，我们无法采用微积分这样的工具来构建最优性条件，因此并没有像连续优化问题那样有比较好的最优性条件(例如 KKT条件)，这也正是造成整数规划问题求解困难的主要原因。基于这个思路，我们尝试从最基本的最优性原理出发来构造整数规划问题的近似最优性条件。

考虑如下的整数规划问题IP：v∗=max⁡{cTx:x∈S⊆Zn}v^*=\max \left\{ c^Tx\ :\ x\in S\subseteq Z^n \right\} v∗=max{cTx : x∈S⊆Zn}
其中 v∗v^{*}v∗为上述整数规划问题的最优解。
假设我们依据某个算法可以得到如下关于该整数规划问题上界的序列：
vˉ1≥vˉ2≥vˉ3≥...≥v∗\bar{v}_1\geq\bar{v}_2\geq\bar{v}_3\geq...\geq v^*vˉ1≥vˉ2≥vˉ3≥...≥v∗
同时我们还可以得到关于该整数规划问题下界的序列：
v‾1≤v‾2≤v‾3≤...≤v∗\underline{v}_{1}\leq\underline{v}_{2}\leq\underline{v}_{3}\leq...\leq{v}^{*} v1≤v2≤v3≤...≤v∗
若对于某个kkk 使得如下表达式成立：vˉk−v‾k≤ϵ\bar{v}_{k}-\underline{v}_{k}\leq\epsilonvˉk−vk≤ϵ，则可得：
vˉk−ϵ≤v∗≤vˉk\bar{v}_{k}-\epsilon \leq v^*\leq\bar{v}_{k}vˉk−ϵ≤v∗≤vˉk

定理1.1 ：若 ϵ=0\epsilon = 0ϵ=0 ，则表明 vˉk−v‾k=0\bar{v}_{k}-\underline{v}_{k} = 0vˉk−vk=0，进一步可得整数规划问题P的最优解为 vˉk\bar{v}_{k}vˉk。

由此可得可行解 vˉk\bar{v}_{k}vˉk 离最优解的距离不超过 ϵ\epsilonϵ。从上面的推导可以知道，如果我们有2个算法分别能够求出原整数规划问题的上界和下界，然后逐步逼近上下界，当上下界之间的距离足够接近的时候，我们就可以得到一个近似最优解。

从上面的分析我们可以看出我们需要分别寻找整数规划问题IP的上界和下界，寻找上界自然我们首先会相当对原整数规划问题做松弛，如下所示：

考虑整数规划问题P的线性松弛问题LP：
u∗=max⁡{cTx:x∈S⊆Rn}u^*=\max \left\{ c^Tx\ :\ x\in S\subseteq R^n \right\} u∗=max{cTx : x∈S⊆Rn}
寻找下界，我们一般会选择对原整数规划问题进行分解：
设原整数规划问题的可行域可以被分解为如下若干个集合的并 S=S1∪S2⋯∪SMS=S_1\cup S_2\cdots \cup S_M S=S1∪S2⋯∪SM我们在每个SiS_iSi小集合上有，zk=max⁡{cTx:x∈Sk⊆Zn}z^k=\max \left\{ c^Tx\ :\ x\in S_k\subseteq Z^n \right\} zk=max{cTx : x∈Sk⊆Zn}
对于k=1,...,Mk=1,...,Mk=1,...,M进一步可得v∗=max⁡kzk≥zkv^* =\underset{k}{\max}z^k\ge z^k v∗=kmaxzk≥zk
将上面两条结合起来松弛问题可以帮助我们得到上界，分解问题可以帮助对整数规划问题进行拆分，同时也可以帮助我们得到下界。把这两个思想结合在一起就构成了我们今天的主角：分支定界法。

2 分支定界法如何实现节省计算量

以如下三个决策变量的整数规划问题为例：
max⁡x1+x2+x3s.t.x1,x2,x3∈{0,1}\max x_1+x_2+x_3\\ s.t.\ x_1,x_2,x_3\in \left\{ 0,1 \right\} maxx1+x2+x3s.t. x1,x2,x3∈{0,1}
上述问题总的可行域为 S={0,1}3S=\left\{0,1 \right\}^3S={0,1}3我们将集合SSS分为两部分S0={x∈S:x1=0}S_0=\left\{ x \in S:x_1=0\right\}S0={x∈S:x1=0}和 S1={x∈S:x1=1}S_1=\left\{ x \in S:x_1=1\right\}S1={x∈S:x1=1}，依次类推S00,S01,...S_{00},S_{01},...S00,S01,...等，如下图所示可以构成一个树结构：

固然我们可以去遍历整个树上的每个节点就可以完成对整数规划问题的求解，但是这样的方法实际上是不可行的，因为其计算时间会随着问题规模急剧增大，这一点在前面我们已经讲过了。所以我们需要通过定出问题的上下界来比较“智能”的给出哪些节点肯定不是最优解，那么就不用搜索这些节点了以便于节省时间。

首先我们需要如下两个定理使得我们的分支定界法变得比较“智能”

定理2.1：（1）若LP不可行，则IP也不可行。（2）设 x∗x^*x∗是LP的最优解，同时有x∗∈Znx^*\in Z^nx∗∈Zn，则x∗x^*x∗也是IP的最优解。
证明：(1)采用反证法直接可以证明。(2) 由于x∗x^*x∗既是IP问题的下界，又是IP问题的上界，所以其必为IP问题的最优解。

定理2.2：设原整数规划问题的可行域可以被分解 S=S1∪S2⋯∪SMS=S_1\cup S_2\cdots \cup S_M S=S1∪S2⋯∪SM 我们在每个 SiS_iSi小集合上有， zk=max⁡{cTx:x∈Sk⊆Zn}z^k=\max \left\{ c^Tx\ :\ x\in S_k\subseteq Z^n \right\} zk=max{cTx : x∈Sk⊆Zn} 对于 k=1,...,Mk=1,...,Mk=1,...,M，设 zˉk\bar{z}^kzˉk为zkz^kzk的一个上界，z‾k\underline{z}^kzk为zkz^kzk的一个下界。由此可得，zˉ=max⁡kzˉk\bar{z}=\underset{k}{\max}\bar{z}^kzˉ=kmaxzˉk是原整数规划问题最优解v∗v^*v∗的上界， z‾=max⁡kz‾k\underline{z}=\underset{k}{\max}\underline{z}^kz=kmaxzk是原整数规划问题的下界。

情况1：Pruned by optimality

如上图所示，图中节点旁边的两个数字分别表示该节点处的上界和下界。根据定理2.2可以对zˉ\bar{z}zˉ进行更新，即 zˉ=max⁡kzˉk=max⁡{20,25}=25\bar{z}=\underset{k}{\max}\bar{z}^k=\max\left\{ 20,25 \right\}=25zˉ=kmaxzˉk=max{20,25}=25，同时也可以对z‾\underline{z}z进行更新，即 z‾=max⁡kz‾k=max⁡{15,20}=20\underline{z}=\underset{k}{\max}\underline{z}^k = \max\left\{ 15,20 \right\}=20z=kmaxzk=max{15,20}=20。我们观察到节点S1S_1S1 处的上界和下界都是20，根据定理1.1可知此时已经求解到了S1S_1S1 问题的最优解了，那么就无需对S1S_1S1 的子节点进一步搜索了。

情况2：Pruned by bound

如上图所示，同样的根据定理2.2可以对zˉ\bar{z}zˉ进行更新，即 zˉ=max⁡kzˉk=max⁡{20,26}=26\bar{z}=\underset{k}{\max}\bar{z}^k=\max\left\{ 20,26\right\}=26zˉ=kmaxzˉk=max{20,26}=26,同时也可以对z‾\underline{z}z
进行更新，即z‾=max⁡kz‾k=max⁡{18,21}=21\underline{z}=\underset{k}{\max}\underline{z}^k = \max\left\{ 18,21 \right\}=21z=kmaxzk=max{18,21}=21。我们观察节点 S1S_1S1可知S1S_1S1的上界是20，而z‾=21\underline{z}=21z=21，也就是说最优解一定是大于等于21的，而S1S_1S1节点的子问题最大只能取到20，由此可知最优解一定不在S1S_1S1处，那么就无需对 S1S_1S1的子节点进一步搜索了。

情况3：No pruning possible

如上图所示，同样的根据定理2.2可以对zˉ\bar{z}zˉ进行更新，zˉ=max⁡kzˉk=max⁡{24,37}=37\bar{z}=\underset{k}{\max}\bar{z}^k=\max\left\{ 24,37\right\}=37zˉ=kmaxzˉk=max{24,37}=37，同时也可以对z‾\underline{z}z进行更新，即z‾=max⁡kz‾k=max⁡{13,−∞}=13\underline{z}=\underset{k}{\max}\underline{z}^k = \max\left\{ 13,-\infty \right\}=13z=kmaxzk=max{13,−∞}=13，观察S1S_1S1和S2S_2S2我们发现不能通过情况1和情况2来跳过节点，也就是说此时我们还得进一步去搜索S1S_1S1和S2S_2S2。

情况 4： Pruned by infeasibility
infeasibility 的情况最简单，就是解到某个节点发现该节点的问题是一个不可行问题，那么必然这个节点之下的子问题也都是不可行的，那么该节点就不用在搜索下去了。

至此，我们总结一下在树搜索过程中会有四种情况发生：

情况1：Pruned by optimality，这种情况是非常幸运的一下子就解到的子问题的最优解。

情况2：Pruned by bound，这种情况可以通过最优解的上下界来排除掉一部分不可能是最优解的节点。

情况3：No pruning possible，这种情况实际上就是最差的情况。

情况4：Pruned by infeasibility ，这种情况因为不可行的原因同样可以删除掉一部分节点不用搜索。了解了分支定界法的四种情况，接下来我们给出分支定界法的算法总流程。

3 分支定界法流程

分支定界算法的流程如下所示：

从图中可以看到初始化阶段我们需要给定输出的全局的上界zˉ\bar{z}zˉ和下界z‾\underline{z}z，如果能有一些启发式的方法获得稍微好点的上下界作为初始解导入那是最好的不过的了。如果没有的话可以先设置为正负无穷大。

接着进入到主循环中，我们通过求解整数规划的连续松弛问题（线性规划）来得到该子问题的上界。在主循环中，剩余的操作其实就是依据我们上一节所讲的四种情况来分别进行，即 Pruned by optimality，Pruned by bound，No pruning possible, Pruned by infeasibility，对于情况1,2,4我们需要提前跳出循环完成剪枝的过程。如果情况1,2,4都不满足的话，就会一直运行到最后一步 return two subproblems 因为没有剪枝，所以需要把两个子节点都添加进来。

4 分支定界法求解0-1线性整数规划问题代码实现

依据上述算法流程，我们在Python中实现了一个基本的分支定界算法。注意的是我们的代码针对是0-1线性整数规划问题，并且目标函数是极大化的情况。同时整数规划模型的建立要依赖Gurobi，所以想要运行代码需要安装Gurobi 以便于代码结果展示。
接下来我们对代码的主干部分做一个简要介绍，首先我们通过Gurobi生成我们要求解的整数规划问题的模型作为我们分支定界法的输入。这部分是Gurobi的指令（不熟悉的Gurobi的同学需要先了解一下Gurobi），如下所示

model = gp.Model("mip1")
x = model.addVar(vtype=GRB.BINARY, name="x")
y = model.addVar(vtype=GRB.BINARY, name="y")
z = model.addVar(vtype=GRB.BINARY, name="z")model.setObjective(x + y + 2 * z, GRB.MAXIMIZE)
model.addConstr(x + 2 * y + 3 * z <= 4, "c0")
model.addConstr(x + y >= 1, "c1")
model.addConstr(x + z == 2, "c2")
model.update()
model.write("model_integer.lp")

准备完输入之后，如下代码实际上就对应我们上一节所述的分支定界算法的主干流程。可以看到 Pruned by optimality，Pruned by bound, Pruned by infeasibility 三种情况会跳出本次循环，而如果三种情况都不满足则就是情况3：No pruning possible，这时候会在本次循环的最后生成两个子问题（child_node1,child_node2）

upper_bound, lower_bound = float('inf'), 0
model_relax = model.relax()
root_node = Node(model = model_relax, upper_bound = upper_bound, lower_bound = lower_bound, candidate_vars = [i for i in range(model.NumVars)])
candidate_node = [root_node]
current_optimum = Nonewhile candidate_node:node, candidate_node = choice_node(candidate_node)if node.upper_bound <= lower_bound:print("prune by bound")continuemodel_status = node.optimize(heuristic_solve)if model_status == 'infeasible':print("prune by infeasiblity")continuenode.update_upper_bound()if node.upper_bound <= lower_bound:print("prune by bound")continueif node.is_integer():node.update_lower_bound()if node.lower_bound > lower_bound:lower_bound = node.lower_boundcurrent_optimum = node.solutioncontinueif node.is_child_problem():child_node1, child_node2 = node.get_child_problem()candidate_node.append(child_node1)candidate_node.append(child_node2)

在代码中有两个模块需要额外提一下，heuristic_solve(problem)函数是负责在每个节点求解线性规划问题的，输入是该节点的问题，输出是该问题的最优解和最优目标函数。在我们的代码中就是直接调用的是Gurobi来帮助我们对线性规划问题进行求解，你当然也可以用其它的方法来求解这个线性规划问题。choice_node(condidate_node)函数是负责每次选择一个节点开始下次计算的。事实上选择节点的方式有非常多的方法，我这里采用的是一个队列的方式来选择节点，即先进先出的规则。节点选择实际上是一个比较大的学问，大家有兴趣的话可以尝试别的选择节点的方式。

def heuristic_solve(problem):problem.Params.OutputFlag = 0problem.optimize()if problem.status == GRB.INFEASIBLE:return None, Nonereturn problem.ObjVal, problem.getVars()def choice_node(condidate_node):node = condidate_node.pop(0)return node, condidate_node

分支定界算法的完整代码见：
https://github.com/WenYuZhi/EasyIntegerProgramming/tree/master/BranchBound
如果大家在使用分支定界算法过程中有一些问题或者发现Bug的话欢迎留言联系我。

5 总结

我们测试了两个简单的整数规划问题，发现我们算法的求解结果和调用Gurobi求解结果是一致的，但是需要指出的是如果决策变量维数稍微大一点的话实际上我们写的分支定界算法也是很慢很慢的，而且随着遍历节点的越来越多，占据的内存空间也非常巨大。
实际上，一方面分支定界作为最最最基本的一种求解整数规划的方法，其效果属实是一般般。真正在工程实际中比较常用的是branch and cut，branch and price，列生成等等“高级”算法。另外一方面，就代码实现角度来说我们写自己写的很土的分支定界算法在一般的问题上肯定是很难超过Gurobi的。我们要知道我们自己撸代码的目的更多是为了理解算法思想，以及为下一步高级算法的设计打下基础。

参考文献：

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