整数规划：分支定界法

胡拉哥

已于 2025-04-15 21:03:59 修改

阅读量7.7k

点赞数 13

分类专栏：优化算法文章标签：算法整数规划运筹学分支定界 python

于 2021-09-03 22:20:08 首次发布

本文链接：https://blog.csdn.net/qx3501332/article/details/120091174

版权

优化算法专栏收录该内容

21 篇文章

订阅专栏

求解整数规划的分支定界算法是一个树搜索（Tree-Search）算法。它的基本思想是考虑整数规划问题的松弛问题，通过增加“整性”约束条件来构造子问题，求解所有子问题从而得到最优的整数解。通过“定界”，可以避免不必要的搜索。

分支定界法

我们用一个例子来说明。考虑下面的整数规划问题：
$\begin{aligned} \max~ & 4x_1+5x_2\\ \text{s.t. } & x_1+4x_2 \leq 10\\ & 3x_1-4x_2 \leq 6\\ & x_1, x_2 \in \mathbb{z}^+ \end{aligned}\quad \text{(IP)}$
注意： $x_1, x_2$ 是正整数。如果不考虑整性约束，那么它就是线性规划问题，称为整数规划问题的松弛问题。

分支

考虑松弛问题:
$\begin{aligned} \max~ & 4x_1+5x_2\\ \text{s.t. } & x_1+4x_2 \leq 10\\ & 3x_1-4x_2 \leq 6\\ & x_1, x_2 \geq 0 \end{aligned}\quad \text{(LP0)}$
我们发现最优解 $x^* = (4, 1.5)$ 。此时， $x_2$ 不是整数。接下来考虑两个子问题 $\text{LP}(1)$ 和 $\text{LP}(2)$ ，分别在父问题的基础上增加约束 $x_2\leq 1$ 和 $x_2 \geq 2$ 。
$\begin{aligned} \max~ & 4x_1+5x_2\\ \text{s.t. } & x_1+4x_2 \leq 10\\ & 3x_1-4x_2 \leq 6\\ & \mathbf{x_2 \leq 1}\\ & x_1, x_2 \geq 0 \end{aligned}\quad (\text{LP1})$

$\begin{aligned} \max~ & 4x_1+5x_2\\ \text{s.t. } & x_1+4x_2 \leq 10\\ & 3x_1-4x_2 \leq 6\\ & \mathbf{x_2 \geq 2}\\ & x_1, x_2 \geq 0 \end{aligned}\quad (\text{LP2})$

再分别求解 $\text{LP}1$ 和 $\text{LP2}$ 。重复这样的步骤，直到无可行解或枚举完所有情形（过程如下图所示）。

在这里插入图片描述

定界

上图搜索树中的叶子结点，要么得到整数可行解，要么无可行解。在实际中，不一定要搜索所有情况，因为可以通过判定最优值的下界（最大化问题），从而过滤掉不必要的分支。

考虑某个结点 $i$ ，假设其最优解是整数解。其目标函数值 $\text{obj}_i$ 是原问题 $\text{IP}$ 的下界，即 $\text{OPT} \geq {\text{obj}_i}$ 。考虑某个节点 $j$ ，假设其最优目标函数值 $\text{obj}_j < \text{obj}_i$ ，那么节点 $j$ 及其子节点可以跳过，因为松弛问题的最优目标函数值是对应整数规划的最优目标函数值的上界。

算法实现

我们用Python实现。线性规划求解器采用谷歌的ORTools。

先定义问题的输入： $c\in\mathbb{R}^n$ , $A\in\mathbb{R}^{m\times n}$ , $b\in\mathbb{R}^m$ 是整数规划问题的参数。此外 $\text{lb}, \text{ub}\in \mathbb{R}^n$ 是决策变量的下界和上界。

import numpy as np


class BranchBound(object):
    """ 分支定界法
    注意：
    1、最大化问题。
    2、整数规划。
    3、决策变量非负。
    """

    def __init__(self, c, A, b, lb=None, ub=None):
        """
        :param c: n * 1 vector
        :param A: m * n matrix
        :param b: m * 1 vector
        :param lb: list, lower bounds of x, e.g. [0, 0, 1, ...]
        :param ub: list, upper bounds of x, e.g. [None, None, ...], None 代表正无穷
        """
        self._c = c
        self._A = A
        self._b = b
        self._lb = lb
        self._ub = ub
        self._m = len(A)
        self._n = len(c)
        self._sol = None  # 决策变量的值
        self._obj_val = -np.inf  # 目标函数值

先直接看分支定界算法的实现，采用深度优先的方式进行搜索。

class BranchBound(object):  
  	
    # ...
  	# 其它函数省略
  
    def _solve_dfs(self, v: Node):
        """ 给定节点v，采用深度优先，搜索v的子树。
        """
        status = v.solve()  # 求解节点v对应的线性规划问题
        # 如果没有可行解，则返回
        if status == pywraplp.Solver.INFEASIBLE:
            return
        # 如果有可行解，接下来搜索v的子树
        sol = v.get_sol()  
        val = v.get_obj_val()
        # 判断v的最优目标函数值是否小于下界，如果是则返回（剪枝）
        # self._obj_val 代表当前的下界，初始值为负无穷
        if val < self._obj_val:
            return
        # 否则找到一个非整数的决策变量，进行分支
        branch_ind = self._get_branch_ind(sol)
        if branch_ind is not None:
            # 构造两个新的节点
            c1, c2 = _get_children(v, branch_ind)
            # 递归搜索c1和c2对应的子树
            self._solve_dfs(c1)
            self._solve_dfs(c2)
        # 更新下界
        elif val > self._obj_val:
            self._sol = sol
            self._obj_val = val

然后是每个节点的定义和方法（包括调用ORTools求解线性规划）。

from ortools.linear_solver import pywraplp


class Node(object):
    
    """ LP relaxation model.
    """

    def __init__(self, c, A, b, lb, ub):
        self._c = c
        self._A = A
        self._b = b
        self.lb = lb
        self.ub = ub
        self._sol = None
        self._obj_val = None

    def _create_model(self):
        """ 根据输入，初始化模型。
        """
        model = pywraplp.Solver.CreateSolver('GLOP')
        m, n = len(self._b), len(self._c)
        # 决策变量
        x = [model.NumVar(self.lb[j], self.ub[j], 'x[%d]' % j) for j in range(n)]
        # 约束: Ax <= b
        for i in range(m):
            ct = model.Constraint(-model.infinity(), self._b[i])
            for j in range(n):
                ct.SetCoefficient(x[j], self._A[i][j])
        # 目标
        obj = model.Objective()
        for j in range(n):
            obj.SetCoefficient(x[j], self._c[j])
        obj.SetMaximization()
        return model, x, obj

    def solve(self):
        """ 求解模型
        """
        model, x, obj = self._create_model()
        status = model.Solve()
        if status == model.OPTIMAL or status == model.FEASIBLE:
            self._sol = [x[j].solution_value() for j in range(len(x))]
            self._obj_val = obj.Value()
        return status

    def get_sol(self):
        return self._sol

    def get_obj_val(self):
        return self._obj_val

最后是生成新节点的方法。

from copy import deepcopy


def _get_children(v: Node, branch_ind: int):
    """ 计算分支节点(其中 i = brand_id)：
        x[i] <= d,
        x[i] >= d+1
        只需要重新设置lb[i]和ub[i]
    """
    d = int(v.get_sol()[branch_ind])
    # left child
    # Python对象是指针引用，所以用深拷贝
    left_child = deepcopy(v)
    left_child.ub[branch_ind] = d
    # right child
    right_child = deepcopy(v)
    right_child.lb[branch_ind] = d+1
    return left_child, right_child