整数规划：分支定界法

入门教程（含代码）

胡拉哥

7166人浏览 · 2021-09-03 22:20:08

胡拉哥 · 2021-09-03 22:20:08 发布

求解整数规划的分支定界算法是一个树搜索（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$ 及其子节点可以跳过，因为松弛问题的最优目标函数值是对应整数规划的最优目标函数值的上界。

如下图所示，我们修改一下 $x_2\geq 2$ 对应的节点，假设它的最优目标函数值是17。

这意味着这个分支所有子节点，无论其最优解是否整数，目标函数值都不会超过17。已知一个整数可行解 $x_1, x_2) = (2, 2)$ ，对应的目标函数值是 18，因此我们不需要再考虑 $x_2\geq 2$ 这个节点及其子节点。

算法实现

我们用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