线性规划技巧: Benders Decomposition

胡拉哥

已于 2023-03-20 17:16:21 修改

阅读量1.1w

点赞数 39

分类专栏：优化算法文章标签：算法线性规划对偶问题

于 2020-03-19 22:41:18 首次发布

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

版权

优化算法专栏收录该内容

21 篇文章 55 订阅

订阅专栏

本文介绍的行生成方法，也称为 Benders 分解，由Jacques F. Benders在1962年提出¹。本文介绍的内容基于 Georrion和Graves².

我们先看一个问题，然后从这个问题出发，来介绍行生成的思想。

设施选址问题

某物流公司要开设仓库，从而服务周边的客户。仓库有“开仓成本”。仓库与客户之间，存在“连接成本”。

已知 $m$ 个"候选仓库"和 $n$ 个客户。要决定开哪些仓库，从而覆盖所有客户，目标是使得总成本最小。

数学规划

根据问题描述，先定义参数和决策变量。

再描述目标和约束，从而得到它的数学规划。

注意， $x_{i,j}\in\{0,1\}$ 可以放松成 $x_{i,j}\geq 0$ ，这不会改变原问题的最优解（证明略）。

行生成

如果问题规模不大，可以试着直接求解。如果解不出来，可以把它分解成小问题，然后用迭代的方式，去逼近最优解。

下面介绍如何分解。

注意决策变量，其中 $x_{i,j}$ 是连续变量， $y_i$ 是离散变量。如果固定 $y_i$ ，那么它变成连续问题，如下所示。

其中 $P_y(x)$ 是线性规划，解起来相对容易。如果把它作为子问题，单独求解行不行？

那么原问题，可以写成这样。

但是，这么做有个麻烦。

子问题的可行域，依赖主问题的 $y$ 值。而主问题的可行域，也依赖子问题的 $x$ 值。两者互相依赖，子问题就不好拆。

下面换个思路。

考虑 $P_y(x)$ 的对偶问题（怎么写对偶？可以看这篇《如何写对偶问题》）。

注意 $y$ 的位置，从约束变成了目标。换句话说，对偶问题的可行域，不受 $y$ 值的影响。

再看原问题，它变成了这样。

引入中间变量 $z$ ，写成如下形式。

看这个形式，似乎更复杂了。约束数量等于 $∣ I ∣$ 。换句话说，约束有无数多个。这更难计算了。

但是，它有个好处。

哪怕只考虑一部分约束，它的解始终是原问题的可行解（因为 $\alpha,\beta$ 可行）。

换句话说，计算的时候，可以只考虑一部分约束，先求得可行解。然后用迭代的方式，逐步逼近最优解。

这样一来，主问题的约束，可以逐个增加。通过求解子问题 $D_y(\alpha,\beta)$ ，得到可行的 $\alpha,\beta$ 值，也是主问题的约束。

每迭代一轮，主问题增加一行，这就是 Benders 分解，也称为行生成方法。

迭代过程

下面介绍具体的迭代过程。

第一步，初始化，令 $z = 0$ ， $I=\emptyset$ 。

主问题的形式如下：

第二步，求解主问题，得到 $y^*, z^*$ 。

注意到，主问题考虑的约束，是原问题约束的子集。因此，其最优目标函数值，是“原问题”的最优目标函数值 $\text{OPT}$ 的下界，即

第三步，求解子问题 $D_{y^*}(\alpha,\beta)$ ，得到 $\alpha^*, \beta^*$ 。

接着，主问题新增一行，即 $I=I\cup \{(\alpha^*, \beta^*)\}$ 。

令 $\text{OPT}_S$ 为子问题的最优目标函数值。回顾原问题的等价形式 $P_1(y,\alpha,\beta)$ 和 $P_2(y, z)$ ，我们有

从而得到 $\text{OPT}$ 的上界，即

第四步，判断停止条件。给定 $\epsilon >0$ ，比如 $\epsilon=1e-6$ ，当 $\text{ub}-\text{lb} < \epsilon$ 时，停止；否则执行第二步。

以上就是行生成方法的迭代过程。

但是，没有完。

按上述步骤，其实解不出来。会出现子问题无解，从而无法迭代。所以需要打个补丁，下面讲怎么打。

可行性

要顺利迭代，就要保证子问题 $D(\alpha,\beta)$ 有可行解。而 $D(\alpha,\beta)$ 是 $P_y(x)$ 的对偶问题，所以 $P_y(x)$ 要有可行解。

换句话说，给定任意的 $y$ ，存在 $x$ 使得下面的约束成立。

容易发现，当 $y = 0$ 时 $x = 0$ ，等式不能满足，此时 $P_y(x)$ 无可行解。

我们理解一下，当 $y = 0$ 时，意味着不开仓。而上面的等式，要求所有客户被连接。不开仓就连接不了客户，所以不可行。

为了连接所有客户，至少要开一个仓。即
$\sum_{i=1}^m y_i \geq 1.$

那么主问题，应该是这样。

这样一来，它解出来的 $y$ 值，使得子问题可行。也就是说，上面的迭代过程，能顺利执行。

最后，还有个小问题。

主问题解的是 $y, z$ ，子问题解的是 $\alpha,\beta$ ，而原问题解的是 $x, y$ 。

已知主问题和子问题的解，怎么得到原问题的解？

回顾子问题 $D_y(\alpha,\beta)$ ，它的对偶问题是 $P_y(x)$ 。换句话说，已知 $\alpha,\beta$ ，它的对偶变量就是 $x$ 。

利用求解器（比如OR-tools），可直接返回 $\alpha,\beta$ 的对偶变量值 $x$ 。

总结

列生成方法，适合求解这种形式的规划：

其中与 $y$ 相关的部分，可以是非线性的。
给定 $y$ ，关于 $x$ 的问题，就是线性规划。考虑它的对偶，从而把与 $y$ 相关的部分，从约束转移到目标。这样一来，原问题的约束，来自对偶问题的解。

那么原问题就可以拆成两个问题：一个是原问题的目标，作为主问题，另一个是原问题的约束，来自对偶问题，作为子问题。

交替求解主问题和子问题，从而逼近原问题的最优解。

值得注意的是，要保证迭代是可行的，就要保证子问题可行。一般来说，只要问题本身合理，可行性就容易满足。

但是可能存在，主问题的某些 $y$ 值导致子问题不可行。添加相应的约束，可以避免这种情况的发生。

Python实现

主问题模型

# master_model.py

from ortools.linear_solver import pywraplp

import numpy as np


class MasterModel(object):

    def __init__(self, f, fix_z=None):
        """
        :param f: 开仓成本(m维向量)
        :param fix_z: 限定z的值. 目的是初始化的时候令z=0,然后求解主问题得到y的值.
        """
        self._solver = pywraplp.Solver('MasterModel',
                                       pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
        self._f = f
        self._fix_z = fix_z
        self._m = len(self._f)
        self._y = None  # 决策变量
        self._z = None  # 决策变量
        self._solution_y = None  # 计算结果
        self._solution_z = None  # 计算结果
        # 迭代求解之前会调用add_constraint
        # 所以决策变量的初始化要放在前面
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()

    def _init_decision_variables(self):
        self._y = [self._solver.NumVar(0, 1, 'y[%d]' % i)
                   for i in range(self._m)]
        self._z = self._solver.NumVar(-self._solver.Infinity(), self._solver.Infinity(), 'z')

    def _init_objective(self):
        obj = self._solver.Objective()
        for i in range(self._m):
            obj.SetCoefficient(self._y[i], self._f[i])
        obj.SetCoefficient(self._z, 1)
        obj.SetMinimization()

    def add_constraint(self, alpha, beta):
        """ Benders主流程会调用此方法为主问题增加新的约束.

        :param alpha: 子问题的解
        :param beta: 子问题的解
        """
        ct = self._solver.Constraint(sum(alpha), self._solver.Infinity())
        for i in range(self._m):
            ct.SetCoefficient(self._y[i], sum(beta[i]))
        ct.SetCoefficient(self._z, 1)

    def _init_constraints(self):
        if self._fix_z is not None:
            ct = self._solver.Constraint(self._fix_z, self._fix_z)
            ct.SetCoefficient(self._z, 1)

        ct = self._solver.Constraint(1, self._solver.infinity())
        for i in range(self._m):
            ct.SetCoefficient(self._y[i], 1)

    def solve(self):
        self._solver.Solve()
        self._solution_y = [self._y[i].solution_value() for i in range(self._m)]
        self._solution_z = self._z.solution_value()

    def get_solution(self):
        return self._solution_y, self._solution_z

    def get_objective_value(self):
        return sum(np.array(self._solution_y) * np.array(self._f)) + self._solution_z

子问题模型

# sub_model.py

from ortools.linear_solver import pywraplp
import numpy as np


class SubModel(object):

    def __init__(self, y, C):
        """
        :param y: 主问题的解(代表y_i=1代表开设仓i)
        :param C: 连接成本(m*n维矩阵)
        """
        self._solver = pywraplp.Solver('SubModel',
                                       pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
        self._y = y
        self._c = C
        self._m = len(self._c)
        self._n = len(self._c[0])
        self._alpha = None  # 决策变量
        self._beta = None  # 决策变量
        self._constraints = []  # 所有约束(后续获取对偶变量, 得到原始问题的解)
        self._solution_alpha = None  # 计算结果
        self._solution_beta = None  # 计算结果

    def _init_decision_variables(self):
        self._alpha = [self._solver.NumVar(-self._solver.Infinity(),
                                           self._solver.Infinity(),
                                           'alpha[%d]' % j)
                       for j in range(self._n)]

        self._beta = [[self._solver.NumVar(0, self._solver.Infinity(), 'beta[%d][%d]' % (i, j))
                      for j in range(self._n)] for i in range(self._m)]

    def _init_constraints(self):
        for i in range(self._m):
            self._constraints.append([])
            for j in range(self._n):
                ct = self._solver.Constraint(-self._solver.Infinity(), self._c[i][j])
                ct.SetCoefficient(self._alpha[j], 1)
                ct.SetCoefficient(self._beta[i][j], -1)
                self._constraints[i].append(ct)

    def _init_objective(self):
        obj = self._solver.Objective()
        for j in range(self._n):
            obj.SetCoefficient(self._alpha[j], 1)
        for i in range(self._m):
            for j in range(self._n):
                obj.SetCoefficient(self._beta[i][j], -self._y[i])

        obj.SetMaximization()

    def solve(self):
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()
        self._solver.Solve()
        self._solution_alpha = [self._alpha[j].solution_value() for j in range(self._n)]
        self._solution_beta = [[self._beta[i][j].solution_value()
                                for j in range(self._n)]
                               for i in range(self._m)]

    def get_solution(self):
        return self._solution_alpha, self._solution_beta

    def get_objective_value(self):
        return sum(self._solution_alpha) - sum(np.array(self._y) * np.sum(self._solution_beta, axis=1))

    def get_dual_values(self):
        """ 得到原问题的解: x[i][j]
        """
        return [[self._constraints[i][j].dual_value()
                for j in range(self._n)]
                for i in range(self._m)]

Benders分解流程

# benders_proc.py

import numpy as np

from master_model import MasterModel
from sub_model import SubModel


class BendersProc(object):
    """ Benders分解流程
    """

    def __init__(self, f, C, max_iter=10000):
        """
        :param f: 开仓成本(m维向量)
        :param C: 连接成本(m*n矩阵), 其中m是候选设施的数量, n是客户的数量
        :param max_iter: 最大循环次数
        """
        self._f = f
        self._c = C
        self._iter_times = 0
        self._max_iter = max_iter
        self._status = -1  # -1:执行错误; 0:最优解; 1: 达到最大循环次数
        self._ub = np.inf
        self._lb = -np.inf
        self._solution_x = None  # 计算结果
        self._solution_y = None  # 计算结果

    def _stop_criteria_is_satisfied(self):
        """ 根据上下界判断是否停止迭代
        """
        if self._ub - self._lb < 0.0001:
            self._status = 0
            return True
        if self._iter_times >= self._max_iter:
            if self._status == -1:
                self._status = 1
            return True
        return False

    def solve(self):
        # 初始令z=0. 求解主问题得到y
        master0 = MasterModel(self._f, fix_z=0)
        master0.solve()
        y, z = master0.get_solution()
        # 下面的迭代需要重新生成master对象
        # 因为master0中z=0是约束条件
        master = MasterModel(self._f)
        sub = None
        # 迭代过程
        while not self._stop_criteria_is_satisfied():
            # 求解子问题
            sub = SubModel(y, self._c)
            sub.solve()
            # 更新上界
            fy = sum(np.array(self._f) * np.array(y))
            self._ub = min(self._ub, fy + sub.get_objective_value())
            # 生成主问题的约束
            alpha, beta = sub.get_solution()
            master.add_constraint(alpha, beta)
            master.solve()
            # 更新y和下界
            y, z = master.get_solution()
            self._lb = master.get_objective_value()

            print(">>> iter %d: lb = %.4f, ub = %.4f" % (self._iter_times, self._lb, self._ub))
            self._iter_times += 1

        # 保存结果
        self._solution_x = sub.get_dual_values()
        self._solution_y = y

        status_str = {-1: "error", 0: "optimal", 1: "attain max iteration"}
        print(">>> Terminated. Status:", status_str[self._status])

    def print_info(self):
        print("---- Solution ----")
        res = {}
        for i in range(len(self._f)):
            if self._solution_y[i] == 0:
                continue
            connected_cities = [str(j) for j in range(len(self._c[0]))
                                if self._solution_x[i][j] > 0]
            res[i] = ', '.join(connected_cities)

        for f, c in res.items():
            print("open facility: %d, connected cities: %s" % (f, c))

主函数

# main.py

from benders_proc import BendersProc
from data import f, C  # data.py包含了一个实例


if __name__ == '__main__':
    bp = BendersProc(f, C)
    bp.solve()
    bp.print_info()

完整代码

参考文献

J.F. Benders. Partitioning Procedures for Solving Mixed-Variables Programming Problems. Numerische Mathematik, Vol. 4, 1962. ↩︎
A.M. Geoffrion and G.W. Graves. Multicommodity distribution system design by benders decomposition. Management Science 20, no. 5, 822–844, 22, 1974. ↩︎

胡拉哥

关注

39
点赞
踩
170

收藏

觉得还不错? 一键收藏
25
评论
线性规划技巧: Benders Decomposition

Benders分解由Jacques F. Benders在1962年提出1. 它是一种把线性规划问题分解为小规模子问题的技巧. 通过迭代求解主问题和子问题, 从而逼近原问题的最优解. 与列生成(Column Generation)相比, Benders分解是一种行生成(Row Generation)技巧: 主问题的约束来自子问题的解. 本文介绍的内容基于 Georrion和Graves2.问题描...
复制链接

扫一扫

专栏目录