Python整数规划—分枝定界法

分枝定界法可用于解纯整数或混合的整数规划问题。在本世纪六十年代初由 Land Doig 和 Dakin 等人提出的。由于这方法灵活且便于用计算机求解，所以现在它已是解整数规划的重要方法。目前已成功地应用于求解生产进度问题、旅行推销员问题、工厂选址问题、背包问题及分配问题等。

分枝定界法

1.定义

对有约束条件的最优化问题（其可行解为有限数）的所有可行解空间恰当地进行系统搜索，这就是分枝与定界内容。通常，把全部可行解空间反复地分割为越来越小的子集，称为分枝；并且对每个子集内的解集计算一个目标下界（对于最小值问题），这称为定界。在每次分枝后，凡是界限超出已知可行解集目标值的那些子集不再进一步分枝，这样，许多子集可不予考虑，这称剪枝。这就是分枝定界法的主要思路。

设有最大化的整数规划问题 A ，与它相应的线性规划为问题 B ，从解问题 B 开始，若其最优解不符合 A 的整数条件，那么 B 的最优目标函数必是 A 的最优目标函数 $z^{*}$ 的上界，记作 $\overline{z}$ ；而 A 的任意可行解的目标函数值将是 $z^{*}$ 的一个下界 $\underline{z}$ 。分枝定界法就是将 B 的可行域分成子区域的方法。逐步减小 $\overline{z}$ 和增大 $\underline{z}$ ，最终求到 $z^{*}$ 。现用下例来说明：

2.例题解说

例 3 求解下述整数规划：

$\begin{aligned} &\operatorname{Max} \quad z=40 x_{1}+90 x_{2} \\ &\left\{\begin{array}{l} 9 x_{1}+7 x_{2} \leq 56 \\ 7 x_{1}+20 x_{2} \leq 70 \\ x_{1}, x_{2} \geq 0 \ \end{array}\right. \end{aligned}$ 且为整数

（i）先不考虑整数限制，即解相应的线性规划B，得最优解为：

使用python单纯形表得出：

        objctive        |                   355.88                   |solution        |  4.81  |  1.82  |   0    |   0    |

最优解为： $x_{1}=4.81, x_{2}=1.82, z=355.88$

可见它不符合整数条件。这时 $z$ 是问题 $A$ 的最优目标函数值 $z^{*}$ 的上界，记作 $\overline{z}$ 。而 $x_{1}=0,x_{2}=0$ 显然是问题 $A$ 的一个整数可行解，这时 $z=0$ ，是 $z^{*}$ 的一个下界，记作 $\underline{z}$ ，即 $0\leq z^{*}\leq 356$ 。

（ii）因为 $x_{1},x_{2}$ 当前均为非整数，故不满足整数要求，任选一个进行分枝。设选 1 x 进行分枝，把可行集分成 2 个子集：

$x_{1}\leq [4.81]=4,x_{1}\geq [4.81]+1=5$

因为 4 与 5 之间无整数，故这两个子集的整数解必与原可行集合整数解一致。这一步称为分枝。这两个子集的规划及求解如下：

问题 $B_{1}$ ：

$\begin{aligned} &\quad \operatorname{Max} \quad z=40 x_{1}+90 x_{2} \\ &\left\{\begin{array}{l} 9 x_{1}+7 x_{2} \leq 56 \\ 7 x_{1}+20 x_{2} \leq 70 \\ 0 \leq x_{1} \leq 4, x_{2} \geq 0 \end{array}\right. \end{aligned}$

最优解为： $x_{2}=4.0,x_{2}=2.1,z_{1}=349$

问题 $B_{2}$ ：

$\begin{aligned} & \quad \text { Max } \quad z=40 x_{1}+90 x_{2} \\ &\left\{\begin{array}{l} 9 x_{1}+7 x_{2} \leq 56 \\ 7 x_{1}+20 x_{2} \leq 70 \\ x_{1} \geq 5, x_{2} \geq 0 \end{array}\right. \end{aligned}$

最优解为： $x_{2}=5.0,x_{2}=1.57,z_{1}=341.4$

再定界： $0\leq z^{*}\leq 349$

（iii）对问题 B1 再进行分枝得问题 B11和 B12 ，它们的最优解为

$\begin{array}{ll} B_{11}: & x_{1}=4, x_{2}=2, z_{11}=340 \\ B_{12}: & x_{1}=1.43, \mathrm{x}_{2}=3.00, z_{12}=327.14 \end{array}$

再定界： $340\leq z^{*}\leq 341$ ，并将 $B_{12}$ 剪枝。

（iv）对问题 $B_{2}$ 再进行分枝得问题 $B_{21}$ 和 $B_{22}$ ，它们的最优解为

$B_{21}:x_{1}=5.44,x_{2}=1.00,z_{22}=308$ 。

$B_{22}:$ 无可行解。

将 $B_{21},B_{22}$ 剪枝。

于是可以断定原问题的最优解为：

$x_{1}=4,x_{2}=2,z^{*}=340$

3.数学建模过程：

分枝定界法求解整数规划（最大化）问题的步骤为：

开始，将要求解的整数规划问题称为问题 A ，将与它相应的线性规划问题称为问题 B 。

解问题 B 可能得到以下情况之一：

（a） B 没有可行解，这时 A 也没有可行解，则停止．

（b） B 有最优解，并符合问题 A 的整数条件， B 的最优解即为 A 的最优解，则停止。

（c） B 有最优解，但不符合问题 A 的整数条件，记它的目标函数值为 $\overline{z}$ 。

用观察法找问题 A 的一个整数可行解，一般可取 $x_{j}=0,j=1,...,n,$ 试探，求得其目标函数值，并记作 $\underline{z}$ 。以 $z^{*}$ 试探，求得其目标函数值，并记作 $\underline{z}\leq z^{*}\leq \overline{z}$ 进行迭代。

建模过程

第一步：分枝，在 B 的最优解中任选一个不符合整数条件的变量 $x_{j}$ ，其值为 $b_{j}$ ，以[ $b_{j}$ ] 表示小于 $b_{j}$ 的最大整数。构造两个约束条件

$x_{j}\leq [b_{j}]$ 和 $x_{j}\geq [b_{j}]+1$

将这两个约束条件，分别加入问题 B ，求两个后继规划问题 $B_{1}$ 和 $B_{2}$ 。不考虑整数条件求解这两个后继问题。

定界，以每个后继问题为一分枝标明求解的结果，与其它问题的解的结果中，找出最优目标函数值最大者作为新的上界 $\overline{z}$ 。从已符合整数条件的各分支中，找出目标函数值为最大者作为新的下界 $\underline{z}$ ，若无作用 $\underline{z}$ 不变。

第二步：比较与剪枝，各分枝的最优目标函数中若有小于 z 者，则剪掉这枝，即以后不再考虑了。若大于 $\underline{z}$ ，且不符合整数条件，则重复第一步骤。一直到最后得到 $z^{*}$ = $\underline{z}$ 为止。得最优整数解 $x_{j}^{*}$ ， $j=1,...,n$ 。

4.编程实现

使用python实现分枝定界法算法：

from scipy.optimize import linprog
import numpy as np
from math import floor, ceil
import copyclass Node(object):def __init__(self, x_bounds=[], freeze_var_list=[], index=0, upper_or_lower=0):self._x_bounds = x_boundsself._freeze_var_list = freeze_var_listself._index = indexself._upper_or_lower = upper_or_lowerprint("创建节点: {}".format(index))print('')def freeze_var(self, index, val):self._x_bounds[index][0] = valself._x_bounds[index][1] = valself._freeze_var_list.append(index)def set_lp_res(self, res):self._res = ress = ""for l in range(len(self._res['x'])):if l in self._freeze_var_list:s += "[" + str(self._res['x'][l]) + "]"else:s += " " + str(self._res['x'][l])print("x: ", s)def check_integer_val_solved(self, m):return True if m == len(self._freeze_var_list) else Falseclass BbAlgorithm(object):def __init__(self, c, a_ub, b_ub, x_b, integer_val):self.c = cself.a_ub = a_ubself.b_ub = b_ubself.x_b = x_bself._integer_val = integer_valself.best_solution = float('inf')self.best_node = Noneself.nodes = []self.nodes_solution = []def solve_lp(self, cur_x_b):return linprog(self.c, A_ub=self.a_ub, b_ub=self.b_ub, bounds=cur_x_b)def check_fessible(self, res):if res['status'] == 0:return Trueelif res['status'] == 2:return Falseelse:raise ("问题无界")def add_node(self, node):res = self.solve_lp(node._x_bounds)if self.check_fessible(res) and res['fun'] < self.best_solution:node.set_lp_res(res)self.nodes_solution.append(res['fun'])self.nodes.append(node)if node.check_integer_val_solved(len(self._integer_val)):self.best_solution = res['fun']self.best_node = nodeprint("----------------当前解决方案-------------------")print("x: ", node._res['x'])print("z: ", node._res['fun'])print("---------------------------------------------------\n")print("==> 将节点添加到树列表: ", node._index)print("==> 当前节点: ", self.nodes_solution)print("")return Trueelse:print("==> 节点不可行: ", node._index)print("==> 当前节点: ", self.nodes_solution)print("")return Falsedef del_higher_val_node(self, z_s):del_list = []for i in range(len(self.nodes_solution)):if self.nodes_solution[i] >= z_s:del_list.append(i)s = ""for i in del_list:s += " " + str(self.nodes[i]._index)print("删除节点: ", s)self.nodes = list(np.delete(self.nodes, del_list))self.nodes_solution = list(np.delete(self.nodes_solution, del_list))print("当前节点: ", self.nodes_solution)print("")def del_item(self, index):print("删除节点: ", self.nodes[index]._index)self.nodes = list(np.delete(self.nodes, index))self.nodes_solution = list(np.delete(self.nodes_solution, index))print("当前节点: ", self.nodes_solution)print("")def check_bounds(self, temp_x_b, index, u_or_l):if u_or_l == 1:if self.x_b[index][0] is not None and temp_x_b[index][0] is None:return Falseelif self.x_b[index][0] is None and temp_x_b[index][0] is not None:return Trueelif self.x_b[index][0] is not None and temp_x_b[index][0] is not None:return False if(self.x_b[index][0] > temp_x_b[index][0]) else Trueelif u_or_l == 2:if self.x_b[index][1] is not None and temp_x_b[index][1] is None:return Falseelif self.x_b[index][1] is None and temp_x_b[index][1] is not None:return Trueelif self.x_b[index][1] is not None and temp_x_b[index][1] is not None:return False if(self.x_b[index][1] < temp_x_b[index][1]) else Trueelse:print("界限误差")exit()def run(self):print("####################### 开始 B & B #####################\n")node_count = 0node = Node(copy.deepcopy(self.x_b), [], node_count)node_count += 1res = self.solve_lp(self.x_b)lower = floor(res['x'][self._integer_val[0]])upper = lower + 1lower_node = Node(copy.deepcopy(self.x_b), [], node_count, 1)lower_node.freeze_var(self._integer_val[0], lower)self.add_node(lower_node)node_count += 1upper_node = Node(copy.deepcopy(self.x_b), [], node_count, 2)upper_node.freeze_var(self._integer_val[0], upper)self.add_node(upper_node)node_count += 1while len(self.nodes) > 0:index = np.argmin(self.nodes_solution)x_b = self.nodes[index]._x_boundsfreeze_list = self.nodes[index]._freeze_var_listres = self.nodes[index]._resfreeze_var_index = len(freeze_list)lower = floor(res['x'][self._integer_val[freeze_var_index]])upper = lower + 1lower_node = Node(copy.deepcopy(x_b), copy.deepcopy(freeze_list), node_count, 1)lower_node.freeze_var(self._integer_val[freeze_var_index], lower)self.add_node(lower_node)node_count += 1upper_node = Node(copy.deepcopy(x_b), copy.deepcopy(freeze_list), node_count, 2)upper_node.freeze_var(self._integer_val[freeze_var_index], upper)self.add_node(upper_node)node_count += 1self.del_item(index)self.del_higher_val_node(self.best_solution)print("############################################################")print("")print("######################### 最佳解决方案 #######################")print(self.best_node._res)if __name__ == "__main__":integer_val = [0,1]c = [-40, -90]A = [[9, 7], [7, 20]]b = [56,70]x_bounds = [[0, None] for _ in range(len(c))]bb_algorithm = BbAlgorithm(c, A, b, x_bounds, integer_val)bb_algorithm.run()

输出结果：

####################### 开始 B & B #####################创建节点: 0创建节点: 1x:  [4.0] 2.0999999999901706
==> 将节点添加到树列表:  1
==> 当前节点:  [-348.99999999911535]创建节点: 2x:  [5.0] 1.5714285714280996
==> 将节点添加到树列表:  2
==> 当前节点:  [-348.99999999911535, -341.4285714285289]创建节点: 3x:  [4.0][2.0]
----------------当前解决方案-------------------
x:  [4. 2.]
z:  -340.0
---------------------------------------------------==> 将节点添加到树列表:  3
==> 当前节点:  [-348.99999999911535, -341.4285714285289, -340.0]创建节点: 4==> 节点不可行:  4
==> 当前节点:  [-348.99999999911535, -341.4285714285289, -340.0]删除节点:  1
当前节点:  [-341.4285714285289, -340.0]删除节点:   3
当前节点:  [-341.4285714285289]############################################################
创建节点: 5==> 节点不可行:  5
==> 当前节点:  [-341.4285714285289]创建节点: 6==> 节点不可行:  6
==> 当前节点:  [-341.4285714285289]删除节点:  2
当前节点:  []删除节点:
当前节点:  []##################################################################################### 最佳解决方案 #######################con: array([], dtype=float64)fun: -340.0message: 'The solution was determined in presolve as there are no non-trivial constraints.'nit: 0slack: array([6., 2.])status: 0success: Truex: array([4., 2.])