• Python整数规划—分枝定界法

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

    分枝定界法

    1.定义

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

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

    2.例题解说

    例 3 求解下述整数规划：

     Maxz=40x1+90x2{9x1+7x2≤567x1+20x2≤70x1,x2≥0 且为整数 

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

    使用python单纯形表得出：

    1.         objctive        |                   355.88                   |
    2.         solution        |  4.81  |  1.82  |   0    |   0    |

    最优解为：x1=4.81,x2=1.82,z=355.88

                            可见它不符合整数条件。这时 z是问题A的最优目标函数值 z∗ 的上界，记作 ¯z 。而 x1=0,x2=0显然是问题 A 的一个整数可行解，这时z=0 ，是 z∗的一个下界，记作 z_， 即0≤z∗≤356。

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

    x1≤[4.81]=4,x1≥[4.81]+1=5

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

    问题B1：

    Maxz=40x1+90x2{9x1+7x2≤567x1+20x2≤700≤x1≤4,x2≥0

    最优解为：x2=4.0,x2=2.1,z1=349

    问题B2：

     Max z=40x1+90x2{9x1+7x2≤567x1+20x2≤70x1≥5,x2≥0

     最优解为：x2=5.0,x2=1.57,z1=341.4

     再定界：0≤z∗≤349

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

    B11:x1=4,x2=2,z11=340B12:x1=1.43,x2=3.00,z12=327.14

    再定界：340≤z∗≤341，并将B12剪枝。

            （iv）对问题 B2再进行分枝得问题 B21和B22，它们的最优解为

    B21:x1=5.44,x2=1.00,z22=308。

                                            ​​​​​​​  B22:  无可行解。

    将 B21,B22 剪枝。

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

    x1=4,x2=2,z∗=340
     

    3.数学建模过程：

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

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

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

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

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

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

    用观察法找问题 A 的一个整数可行解，一般可取xj=0,j=1,...,n,  试探， 求得其目标函数值，并记作z_。以 z∗​ ​​​​​​试探， 求得其目标函数值，并记作z_≤z∗≤¯z进行迭代。

    建模过程

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

     xj≤[bj]和xj≥[bj]+1

    将这两个约束条件，分别加入问题 B ，求两个后继规划问题 B1 和 B2 。不考虑整数条件 求解这两个后继问题。

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

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

     

     

     4.编程实现

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

    1. from scipy.optimize import linprog
    2. import numpy as np
    3. from math import floor, ceil
    4. import copy
    5.  
    6.  
    7. class Node(object):
    8.     def __init__(self, x_bounds=[], freeze_var_list=[], index=0, upper_or_lower=0):
    9.         self._x_bounds = x_bounds
    10.         self._freeze_var_list = freeze_var_list
    11.         self._index = index
    12.         self._upper_or_lower = upper_or_lower
    13.  
    14.         print("创建节点: {}".format(index))
    15.         print('')
    16.  
    17.     def freeze_var(self, index, val):
    18.         self._x_bounds[index][0] = val
    19.         self._x_bounds[index][1] = val
    20.         self._freeze_var_list.append(index)
    21.  
    22.     def set_lp_res(self, res):
    23.         self._res = res
    24.         s = ""
    25.         for l in range(len(self._res['x'])):
    26.             if l in self._freeze_var_list:
    27.                 s += "[" + str(self._res['x'][l]) + "]"
    28.             else:
    29.                 s += " " + str(self._res['x'][l])
    30.         print("x: ", s)
    31.  
    32.     def check_integer_val_solved(self, m):
    33.         return True if m == len(self._freeze_var_list) else False
    34.  
    35.  
    36. class BbAlgorithm(object):
    37.     def __init__(self, c, a_ub, b_ub, x_b, integer_val):
    38.         self.c = c
    39.         self.a_ub = a_ub
    40.         self.b_ub = b_ub
    41.         self.x_b = x_b
    42.         self._integer_val = integer_val
    43.         self.best_solution = float('inf')
    44.         self.best_node = None
    45.         self.nodes = []
    46.         self.nodes_solution = []
    47.  
    48.     def solve_lp(self, cur_x_b):
    49.         return linprog(self.c, A_ub=self.a_ub, b_ub=self.b_ub, bounds=cur_x_b)
    50.  
    51.     def check_fessible(self, res):
    52.         if res['status'] == 0:
    53.             return True
    54.         elif res['status'] == 2:
    55.             return False
    56.         else:
    57.             raise ("问题无界")
    58.  
    59.     def add_node(self, node):
    60.         res = self.solve_lp(node._x_bounds)
    61.         if self.check_fessible(res) and res['fun'] < self.best_solution:
    62.             node.set_lp_res(res)
    63.             self.nodes_solution.append(res['fun'])
    64.             self.nodes.append(node)
    65.             if node.check_integer_val_solved(len(self._integer_val)):
    66.                 self.best_solution = res['fun']
    67.                 self.best_node = node
    68.                 print("----------------当前解决方案-------------------")
    69.                 print("x: ", node._res['x'])
    70.                 print("z: ", node._res['fun'])
    71.                 print("---------------------------------------------------\n")
    72.             print("==> 将节点添加到树列表: ", node._index)
    73.             print("==> 当前节点: ", self.nodes_solution)
    74.             print("")
    75.             return True
    76.         else:
    77.             print("==> 节点不可行: ", node._index)
    78.             print("==> 当前节点: ", self.nodes_solution)
    79.             print("")
    80.             return False
    81.  
    82.     def del_higher_val_node(self, z_s):
    83.         del_list = []
    84.         for i in range(len(self.nodes_solution)):
    85.             if self.nodes_solution[i] >= z_s:
    86.                 del_list.append(i)
    87.         s = ""
    88.         for i in del_list:
    89.             s += " " + str(self.nodes[i]._index)
    90.         print("删除节点: ", s)
    91.         self.nodes = list(np.delete(self.nodes, del_list))
    92.         self.nodes_solution = list(np.delete(self.nodes_solution, del_list))
    93.         print("当前节点: ", self.nodes_solution)
    94.         print("")
    95.  
    96.     def del_item(self, index):
    97.         print("删除节点: ", self.nodes[index]._index)
    98.         self.nodes = list(np.delete(self.nodes, index))
    99.         self.nodes_solution = list(np.delete(self.nodes_solution, index))
    100.         print("当前节点: ", self.nodes_solution)
    101.         print("")
    102.  
    103.     def check_bounds(self, temp_x_b, index, u_or_l):
    104.         if u_or_l == 1:
    105.             if self.x_b[index][0] is not None and temp_x_b[index][0] is None:
    106.                 return False
    107.             elif self.x_b[index][0] is None and temp_x_b[index][0] is not None:
    108.                 return True
    109.             elif self.x_b[index][0] is not None and temp_x_b[index][0] is not None:
    110.                 return False if(self.x_b[index][0] > temp_x_b[index][0]) else True
    111.         elif u_or_l == 2:
    112.             if self.x_b[index][1] is not None and temp_x_b[index][1] is None:
    113.                 return False
    114.             elif self.x_b[index][1] is None and temp_x_b[index][1] is not None:
    115.                 return True
    116.             elif self.x_b[index][1] is not None and temp_x_b[index][1] is not None:
    117.                 return False if(self.x_b[index][1] < temp_x_b[index][1]) else True
    118.         else:
    119.             print("界限误差")
    120.             exit()
    121.  
    122.     def run(self):
    123.         print("####################### 开始 B & B #####################\n")
    124.         node_count = 0
    125.         node = Node(copy.deepcopy(self.x_b), [], node_count)
    126.         node_count += 1
    127.         res = self.solve_lp(self.x_b)
    128.  
    129.         lower = floor(res['x'][self._integer_val[0]])
    130.         upper = lower + 1
    131.  
    132.         lower_node = Node(copy.deepcopy(self.x_b), [], node_count, 1)
    133.         lower_node.freeze_var(self._integer_val[0], lower)
    134.         self.add_node(lower_node)
    135.         node_count += 1
    136.  
    137.         upper_node = Node(copy.deepcopy(self.x_b), [], node_count, 2)
    138.         upper_node.freeze_var(self._integer_val[0], upper)
    139.         self.add_node(upper_node)
    140.         node_count += 1
    141.  
    142.         while len(self.nodes) > 0:
    143.             index = np.argmin(self.nodes_solution)
    144.             x_b = self.nodes[index]._x_bounds
    145.             freeze_list = self.nodes[index]._freeze_var_list
    146.             res = self.nodes[index]._res
    147.             freeze_var_index = len(freeze_list)
    148.  
    149.             lower = floor(res['x'][self._integer_val[freeze_var_index]])
    150.             upper = lower + 1
    151.             lower_node = Node(copy.deepcopy(x_b), copy.deepcopy(freeze_list), node_count, 1)
    152.             lower_node.freeze_var(self._integer_val[freeze_var_index], lower)
    153.             self.add_node(lower_node)
    154.             node_count += 1
    155.  
    156.             upper_node = Node(copy.deepcopy(x_b), copy.deepcopy(freeze_list), node_count, 2)
    157.             upper_node.freeze_var(self._integer_val[freeze_var_index], upper)
    158.             self.add_node(upper_node)
    159.             node_count += 1
    160.  
    161.             self.del_item(index)
    162.             self.del_higher_val_node(self.best_solution)
    163.             print("############################################################")
    164.         print("")
    165.         print("######################### 最佳解决方案 #######################")
    166.         print(self.best_node._res)
    167.  
    168.  
    169. if __name__ == "__main__":
    170.     integer_val = [0,1]
    171.     c = [-40, -90]
    172.     A = [[9, 7], [7, 20]]
    173.     b = [56,70]
    174.     x_bounds = [[0, None] for _ in range(len(c))]
    175.     bb_algorithm = BbAlgorithm(c, A, b, x_bounds, integer_val)
    176.     bb_algorithm.run()

     输出结果：

    1. ####################### 开始 B & B #####################
    2.  
    3. 创建节点: 0
    4.  
    5. 创建节点: 1
    6.  
    7. x:  [4.0] 2.0999999999901706
    8. ==> 将节点添加到树列表:  1
    9. ==> 当前节点:  [-348.99999999911535]
    10.  
    11. 创建节点: 2
    12.  
    13. x:  [5.0] 1.5714285714280996
    14. ==> 将节点添加到树列表:  2
    15. ==> 当前节点:  [-348.99999999911535, -341.4285714285289]
    16.  
    17. 创建节点: 3
    18.  
    19. x:  [4.0][2.0]
    20. ----------------当前解决方案-------------------
    21. x:  [4. 2.]
    22. z:  -340.0
    23. ---------------------------------------------------
    24.  
    25. ==> 将节点添加到树列表:  3
    26. ==> 当前节点:  [-348.99999999911535, -341.4285714285289, -340.0]
    27.  
    28. 创建节点: 4
    29.  
    30. ==> 节点不可行:  4
    31. ==> 当前节点:  [-348.99999999911535, -341.4285714285289, -340.0]
    32.  
    33. 删除节点:  1
    34. 当前节点:  [-341.4285714285289, -340.0]
    35.  
    36. 删除节点:   3
    37. 当前节点:  [-341.4285714285289]
    38.  
    39. ############################################################
    40. 创建节点: 5
    41.  
    42. ==> 节点不可行:  5
    43. ==> 当前节点:  [-341.4285714285289]
    44.  
    45. 创建节点: 6
    46.  
    47. ==> 节点不可行:  6
    48. ==> 当前节点:  [-341.4285714285289]
    49.  
    50. 删除节点:  2
    51. 当前节点:  []
    52.  
    53. 删除节点:  
    54. 当前节点:  []
    55.  
    56. ############################################################
    57.  
    58. ######################### 最佳解决方案 #######################
    59.      con: array([], dtype=float64)
    60.      fun: -340.0
    61.  message: 'The solution was determined in presolve as there are no non-trivial constraints.'
    62.      nit: 0
    63.    slack: array([6., 2.])
    64.   status: 0
    65.  success: True
    66.        x: array([4., 2.])

  • 相关阅读:
    特殊文件（Properties属性文件）
    C++——类与对象(上)
    python + urllib + BeautifulSoup 获取百度首页标题
    MySQL云数据库5.5导入到自建MySQL数据库5.7
    网络编程_fd_set结构
    2023年云计算的发展趋势如何？
    报名即将结束！11 大云原生领域开源技术干货一场拿下
    windows下python opencv ffmpeg读取摄像头实现rtsp推流 拉流
    C++：详细的说明智能指针的使用以及底层实现，以及删除器和包装器的使用
    【图像融合】基于matlab粒子群优化自适应多光谱图像融合【含Matlab源码 004期】
  • 原文地址：https://blog.csdn.net/qq_21402983/article/details/126388515