分支定界(隐枚举)求解MIP的复现

Word count: 1.4k   |   Reading time≈ 7 min

分支定界求解MIP的复现

算法思想

分支定界法是本科运筹学课程中学习的求解整数规划的算法之一，分支定界算法的核心思想就是分枝和剪枝。当我们不考虑所求解必须是整数这个条件时，用单纯形法可求出最优解，但是这个解往往不全是整数，因此我们采用剪枝的方式一点一点缩小范围，直到所求解为整数解。

  • 分枝定界法是一种基于“树”结构的搜索或遍历方法
  • 分枝定界法中的整数性判断不会在本质上引起数值困难
  • 定界是为了避免无效的分枝搜索，恰当的分枝有助于更好定界
  • 分枝定界法是部分枚举而不是穷举

算法程序结构和代码

算法由初始化节点，更新上界，更新下界，增加约束获得子问题，求解子问题几部分构成。具体就是在算法主循环中维护一个节点队列和一个等待处理的变量队列，每次循环中出队列一个节点，进行初步计算和判定上下界，并进行减枝，如果顺利计算，则从变量队列中出下一个变量进行分枝，将变量中增加cut后的两个节点push入节点队列，进入下一循环。终止条件是需要分枝计算的变量队列为空。

测试代码和算法class如下。具体使用方式如测试代码所示，输入gurobi.model的建立好的模型和模型中的零一量的namelist即可（后面是用getvarbyname处理的，因此输入的是字符串）。

测试代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


import gurobipy as gp
from gurobipy import GRB
from branch_and_bound_class import *


model = gp.Model("mip1")
x = model.addVars(10, name = 'x', vtype = GRB.BINARY)

model.setObjective(-x[0]  -x[1]  -2*x[2]  -2*x[8] - x[9], GRB.MINIMIZE)
model.addConstr(x[0] + 2*x[1] + 3*x[2] + 5*x[3] + 3*x[4] <= 8, "c0")
model.addConstr(2*x[3] + 2*x[4] + 3*x[5] + 5*x[6] + 3*x[7] <= 10, "c1")
model.addConstr(x[7] + x[8] + 3*x[9] <= 4, "c2")
model.addConstr(2*x[0] + x[2] + 3*x[7] + 3*x[8] + 2*x[9] <= 8, "c3")
model.addConstr(x[7] + x[8] + 3*x[9] >= 1, "c4")
model.addConstr(2*x[4] + 2*x[5] + x[6] + 5*x[7] + 3*x[8] >= 4, "c5")
# model.optimize()


candidate_vars = ['x[0]','x[1]','x[2]','x[3]','x[4]','x[5]','x[6]','x[7]','x[8]','x[9]']




solve(model,candidate_vars)

算法class

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
import gurobipy as gp
from gurobipy import GRB

"""branch and bound code
输入gurobi的MIP模型 目前仅仅支持binary variable
"""


class Node:
    def __init__(self, model, upper_bound, lower_bound, candidate_vars):
        """初始化节点

        Args:
            model (_type_): gurobi model object
            upper_bound (_type_): _description_
            lower_bound (_type_): _description_
            candidate_vars (_type_): list,NumVars
        """
        self.upper_bound, self.lower_bound = upper_bound, lower_bound
        self.model = model
        # print(candidate_vars,upper_bound,lower_bound)
        self.candidate_vars = candidate_vars.copy()
        assert(upper_bound >= lower_bound), "upper bound is less than lower bound"

    def optimize(self, heuristic_solve):
        """求解Node

        Args:
            heuristic_solve (anotherFuc): 传入函数参数,求解节点的方法

        Returns:
            _type_: _description_
        """
        self.obj_values, self.solution = heuristic_solve(self.model)
        if self.obj_values == None:
            return "infeasible"
        return "feasible"
    
    def update_upper_bound(self):
        """更新upbound,最小化问题的话,上界无条件更新
        """
        # if self.upper_bound > self.obj_values:
        self.upper_bound = self.obj_values
        assert(self.lower_bound <= self.obj_values)
        assert(self.lower_bound <= self.upper_bound), "upper bound is less than lower bound"

    def update_lower_bound(self):
        """更新lowerbound,最小化问题下界条件更新
        """
        if self.lower_bound < self.obj_values:
            self.lower_bound = self.obj_values
            assert(self.lower_bound <= self.obj_values)
            assert(self.lower_bound <= self.upper_bound), "upper bound is less than lower bound"
        
    def is_integer(self):
        """判断解的质量

        Returns:
            _type_: True则表示解OK
        """
        for var in self.solution:
            if 0 < var.x and var.x < 1:
                return False
        return True
    
    def is_child_problem(self):
        # print(self.candidate_vars)
        # print(len(self.candidate_vars))
        if len(self.candidate_vars) > 0:
            return True
    
    def get_child_problem(self):
        """_summary_得替换成GRBgetvarbyname,因为不是所有变量都是binary

        Returns:
            _type_: _description_
        """
        self.child_left, self.child_right = self.model.copy(), self.model.copy()
        branch_index, self.condidate_child_vars = self.choice_branch(self.candidate_vars)
        # 分枝left bound 和 right bound。
        # print(branch_index)
        # print(self.child_left.getVarByName('x0'))
        # for v in self.child_left.getVars():
        #     print('%s' % (v.VarName,))
        self.child_left.addConstr(self.child_left.getVarByName(branch_index) == 0,"left")
        self.child_right.addConstr(self.child_right.getVarByName(branch_index) == 1,"right")
        self.child_left.write("left.lp")
        self.child_right.write("right.lp")
        node_left = Node(self.child_left, self.upper_bound, self.lower_bound, self.condidate_child_vars)
        node_right = Node(self.child_right, self.upper_bound, self.lower_bound, self.condidate_child_vars)

        return node_left, node_right
    
    def choice_branch(self, candidate_vars):
        """选择分枝的变量,这里可以加一些优化,但是似乎不太好加   现在是从等待分枝的变量中直接pop0

        Args:
            candidate_vars (_type_): 维护一个栈,剩余的需要分枝的binary variables

        Returns:
            _type_: _description_
        """
        self.condidate_child_vars = self.candidate_vars.copy()
        branch_index = self.condidate_child_vars.pop(0)  # 改为name后应该弹出的是str：varName
        return branch_index, self.condidate_child_vars
    
    def write(self):
        self.model.write("model.lp")
    

def heuristic_solve(problem):
    problem.Params.OutputFlag = 0
    problem.optimize()
    if problem.status == GRB.INFEASIBLE:
        return None, None
    # for v in problem.getVars():
    #     print('%s %g' % (v.VarName, v.X))
    # print('Obj: %g' % problem.ObjVal)
    return problem.ObjVal, problem.getVars()

def choice_node(condidate_node):
    """选择下一个计算的节点 现在也是pop0

    Args:
        condidate_node (_type_): 维护的还没计算的节点

    Returns:
        _type_: _description_
    """
    node = condidate_node.pop(0)
    return node, condidate_node


def solve(model,candidate_vars):
    model.update()
    upper_bound, lower_bound = float('inf'), float('-inf')
    model_relax = model.relax()
    root_node = Node(model = model_relax, upper_bound = upper_bound, lower_bound = lower_bound, candidate_vars = candidate_vars)
    candidate_node = [root_node]
    current_optimum = None
    while candidate_node:
        # print(candidate_node)
        node, candidate_node = choice_node(candidate_node)
        if node.upper_bound <= lower_bound:
            print("prune by bound")
            continue
        model_status = node.optimize(heuristic_solve)
        if model_status == 'infeasible':
            print("prune by infeasiblity")
            continue
        node.update_lower_bound()
        if node.lower_bound >= upper_bound:
            print("prune by bound")
            continue
        if node.is_integer():
            # print('yes')
            # exit(0)
            node.update_upper_bound()
            # lower bound 
            if node.upper_bound < upper_bound:
                upper_bound = node.upper_bound
                current_optimum = node.solution
            continue

        if node.is_child_problem():
            child_node1, child_node2 = node.get_child_problem()
            candidate_node.append(child_node1)
            candidate_node.append(child_node2)
    print('lower_bound: ',lower_bound)
    print("upper_bound: ", upper_bound)
    print("optimum:", current_optimum)
Donate
  • Copyright： Copyright is owned by the author. For commercial reprints, please contact the author for authorization. For non-commercial reprints, please indicate the source.
  • Copyrights © 2015-2023 galaxy
  • Visitors: | Views:

请我喝杯咖啡吧~

支付宝
微信