Python实现分枝定界法求解整数规划问题

分枝定界法（Branch and Bound）是一种求解整数规划问题的常用算法，其既可以求解纯整数规划问题，也可以求解混合整数规划问题。这种方法的基本思想是对有约束条件的最优化问题的所有可行解空间进行搜索。

在分枝定界法中，全部可行解空间被反复地分割为越来越小的子集，这个过程称为“分支”。同时，对每个子集内的解集计算一个目标下界（对于最小值问题）或上界（对于最大值问题），这个过程称为“定界”。

在每次分枝后，凡是界限超出已知可行解集目标值的那些子集不再进一步分枝，这样可以节省计算资源，这个过程称为“剪枝”。通过这种方式，算法能够逐步缩小搜索范围，直至找到最优解。

分枝定界法的具体步骤包括：

1. 问题建模：首先，将原问题转化为组合优化问题，并确定优化目标。这是将实际问题抽象化的过程，为后续的算法处理打下基础。
2. 初始化树：将原问题构建成一个树状结构，其中每个节点代表问题的一个子问题。这个树状结构用于表示问题的所有可能解空间，并作为后续分枝和定界的基础。
3. 选择分支：从树中选择一个节点进行分枝。通常，选择一个约束条件，将问题拆分成两个或多个子问题。这些子问题将作为新的节点添加到树中。
4. 求解子问题：对每个子问题，使用适当的方法（如线性规划）求解，得到子问题的解以及对应的目标函数值。
5. 定界：根据求解得到的子问题解，更新整个问题的上界和下界。对于最小化问题，下界是到目前为止找到的最小目标函数值，上界则是可能的最大目标函数值；对于最大化问题则相反。
6. 剪枝：如果某个子问题的目标函数值已经超过当前已知的最优解（对于最小化问题）或者低于当前已知的最劣解（对于最大化问题），那么可以剪去这个子问题的整个树，因为它不可能包含更好的解。剪枝是分枝定界法中减少计算量的关键步骤。
7. 重复过程：继续选择未处理的节点进行分枝、求解子问题、定界和剪枝，直到找到最优解或确定无解为止。

在整个过程中，分枝定界法通过不断地将问题分解为更小的子问题，并利用定界和剪枝技术来减少计算量，从而有效地求解整数规划问题。以下是Python实现分枝定界法：

import gurobipy as gp
from gurobipy import GRB

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

def choice_node(condidate_node):
    node = condidate_node.pop(0)
    return node, condidate_node

class Node:
    def __init__(self, model, upper_bound, lower_bound, candidate_vars):
        self.upper_bound, self.lower_bound = upper_bound, lower_bound
        self.model = model
        self.candidate_vars = candidate_vars.copy()
        assert(upper_bound >= lower_bound), "upper bound is less than lower bound"

    def optimize(self, heuristic_solve):
        self.obj_values, self.solution = heuristic_solve(self.model)
        if self.obj_values == None:
            return "infeasible"
        return "feasible"

    def update_upper_bound(self):
        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):
        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):
        for var in self.solution:
            if 0 < var.x and var.x < 1:
                return False
        return True

    def is_child_problem(self):
        if self.candidate_vars:
            return True

    def get_child_problem(self):
        self.child_left, self.child_right = self.model.copy(), self.model.copy()
        branch_index, self.condidate_child_vars = self.choice_branch(self.candidate_vars)
        self.child_left.addConstr(self.child_left.getVars()[branch_index] == 0)
        self.child_right.addConstr(self.child_right.getVars()[branch_index] == 1)
        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):
        self.condidate_child_vars = self.candidate_vars.copy()
        branch_index = self.condidate_child_vars.pop(0)
        return branch_index, self.condidate_child_vars

    def write(self):
        self.model.write("model.lp")


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.MAXIMIZE)
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()
model.write("model_integer.lp")

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 = None

while 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_upper_bound()
    if node.upper_bound <= lower_bound:
        print("prune by bound")
        continue
    if node.is_integer():
        node.update_lower_bound()
        if node.lower_bound > lower_bound:
            lower_bound = node.lower_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("optimum:", current_optimum)

运行结果如下：