线性规划技巧: Dantzig&Wolfe Decomposti

作者: 胡拉哥 | 来源:发表于2020-03-27 20:48 被阅读0次

Dantzig&Wolfe分解(简称DW分解)^[1]是一种列生成技巧,可以把一类特殊形式线性规划问题分解成若干子问题进行求解.

问题描述

我们考虑如下形式的线性规划问题:

$\begin{aligned} \min ~ & c_1^Tx_1 & + ~ & c_2^Tx_2 & + ~ & \cdots & + ~ & c^T_m x_m & \\ \text{s.t. } & A_{0,1}x_1 & + ~ & A_{0,2}x_2 & + ~ & \cdots & + ~& A_{0,m} x_m & = b_0 \\ & A_{1,1} x_1 & & & & & & & = b_1 \\ & & & A_{2,2} x_2 & & & & & = b_2 \\ & & & & & ... & & & \\ & & & & & & & A_{n,m}x_m & = b_m \\ & x_j \geq 0, & \forall j & \end{aligned}$

其中 $A_{i,j}\in \mathbb{R}^{m_i \times n_j}$ , $b_i \in\mathbb{R}^{m_i}$ , $c_j, x_j \in \mathbb{R}^{n_j}$ .

说明

上述规划共 $m+1$ 组约束. 第一组约束包含了所有的决策变量, 称为连接约束(Linking constraints), 接下来是 $m$ 组独立的约束.
$c_j, x_j, b_i$ 都是向量.
上述规划如果写成标准形式 $\min_x \{c^Tx | Ax = b, x\geq 0\}$ , 其系数矩阵 $A$ 的维度是 $\sum_{i=0}^m m_i \times \sum_{j=1}^n n_j$ . 我们发现 $A$ 实际上是相对稀疏的. 当 $m$ , $n$ 的非常大时, 矩阵的规模可能大到无法直接求解上述规划. 在这种情况下, 我们考虑把它分解成多个子问题进行求解(DW分解).

应用场景

例1 一个公司有 $m$ 个部门, 各部门有独立的约束, 部门之间也有约束. 目标是最小化所有部门的总成本.

例2 一个零售公司有 $m$ 个仓库, 需要决定每个仓库存放的商品. 每个仓库中对商品有一些约束. 商品关于各仓库也需要满足一些约束. 目标是最小化总的成本(例如配送成本/时效等).

主问题

令 $P_i = \{ x_i \in \mathbb{R}^{n_i} \mid A_{i,i}x_i = b_i\}$ . 如果 $P_i$ 是有界的, 那么 $x_i$ 可以表示成 $P_i$ 顶点的凸组合. 设 $v_{i,1}, v_{i,2}, \ldots, v_{i, p_i}$ 代表 $P_i$ 的顶点(Extreme point), 则存在 $\lambda_{i,j} \geq 0$ 且 $\sum_{j=1}^{p_i} \lambda_{i,j}=1$ 使得
$x_i = \sum_{j=1}^{p_i} \lambda_{i,j} v_{i,j}.$
在一般情况下( $P_i$ 无界或有界时), 它可以用顶点和极射线的线性组合来表示(参考 Minkowski表示定理). 具体来说, 令 $w_{i,1}, w_{i,2}, \ldots, w_{i,q_i}$ 代表极射线(Extreme ray). 则存在 $\lambda_{i,j} \geq 0$ 且 $\sum_{j=1}^{p_i} \lambda_{i,j}=1$ , $\mu_{i,j}\geq 0$ 使得
$x_i = \sum_{j=1}^{p_i} \lambda_{i,j} v_{i,j} + \sum_{j=1}^{q_i}\mu_{i,j} w_{i,j}.$

假设我们枚举 $P_1, P_2, \ldots, P_m$ 所有的顶点和极射线. 把上式代入原问题得到主问题(Master problem):

主问题

$\begin{aligned} \min\ & \sum_{i=1}^m\sum_{j=1}^{p_i}\lambda_{i,j}(c_i^Tv_{i,j}) + \sum_{i=1}^m\sum_{j=1}^{q_i}\mu_{i,j}(c_i^Tw_{i,j}) \\ \text{s.t. } & \sum_{i=1}^m \sum_{j=1}^{p_i}\lambda_{i,j} A_{0,i}v_{i,j} + \sum_{i=1}^m \sum_{j=1}^{p_i}\mu_{i,j} A_{0,i}w_{i,j} = b_0 \\ & \sum_{j=1}^{p_i} \lambda_{i,j} =1, \quad \forall i \\ & \lambda_{i,j}, \mu_{i,j} \geq0, \quad \forall i, j. \end{aligned}$

说明

$\lambda_{i,j}$ , $\mu_{i,j}$ 是决策变量.
约束的数量为 $m_0+m$ (第一个等式包含了 $m_0$ 个约束). 原问题的约束数量是 $\sum_{i=0}^m m_i$ . 因此主问题的约束数量明显减少了.
但是主问题的变量显著增加了( $\lambda_{i,j}, \mu_{i,j}$ 对应所有的顶点和极射线). 与列生成技巧相似, 主问题一开始只考虑两个可行解(对应顶点和极射线). 通过求解子问题得到新的顶点或极射线.
在实际问题中, 一般情况下 $P_i$ 都是有界的. 此时主问题可以简化成如下形式.
$\begin{aligned} \min\ & \sum_{i=1}^m\sum_{j=1}^{p_i}\lambda_{i,j}(c_i^Tv_{i,j}) \\ \text{s.t. } & \sum_{i=1}^m \sum_{j=1}^{p_i}\lambda_{i,j} A_{0,i}v_{i,j} = b_0 \\ & \sum_{j=1}^{p_i} \lambda_{i,j} =1, \quad \forall i \\ & \lambda_{i,j} \geq 0, \quad \forall i, j. \end{aligned}$

子问题

定义主问题的对偶变量 $y\in\mathbb{R}^{m_0}$ , $z_i \in \mathbb{R}$ , $i=1, 2, \ldots, m$ . 我们计算变量 $\lambda_{i,j}, \mu_{i,j}$ 的Reduced Cost:

$\begin{aligned} & \alpha_{i,j} = c_i^T v_{i,j} - y^TA_{0,i}v_{i,j} - z_i \\ & \beta_{i, j} = c_i^Tw_{i,j} - y^TA_{0,i}w_{i,j} \end{aligned}$

其中 $\alpha_{i,j}$ , $\beta_{i,j}$ 分别代表 $\lambda_{i,j}, \mu_{i,j}$ 的reduced cost. 当 $\alpha_{i,j}$ , $\beta_{i,j}\geq 0$ 时得到原问题的最优解. 因此在子问题中, 我们需要计算 $v_{i,j}$ (或 $w_{i,j}$ )使得 $\alpha_{i,j}, \beta_{i,j}$ 的值尽可能小.

子问题 - $i$

$\begin{aligned} \min\ & f_{i} := (c_i - A_{0,i}^Ty)^T x_i \\ \text{s.t. } & A_{i,i} x_i = b_i \\ & x_i \geq 0. \end{aligned}$

求解子问题时考虑三种情况.

最优解的目标值为 $-\infty$ . 此时 $\alpha_{i,j}, \beta_{i,j} < 0$ , 子问题的解是极射线 $w_{i,j}$ . 我们把 $A_{0,i}w_{i,j}$ 加入到主问题进行求解.
最优解的目标值有界且 $f_i < z_i$ . 此时 $\alpha_{i,j} < 0$ , 子问题的解是顶点 $v_{i,j}$ . 我们把 $A_{0,i}v_{i,j}$ 加入到主问题进行求解.
最优解的目标值有界且 $f_{i,j} \geq z_i$ . 此时 $0\leq \alpha_{i,j} \leq \beta_{i,j}$ (注意到实际上 $z_i\geq 0$ ), 得到最优解.

例子: 一个选品问题

考虑 $m$ 个零售品牌商, 每个零售品牌商有自己的商品(SKU), 例如可口可乐, 雪碧和芬达对应同一家品牌商. 一家电商平台需要从 $m$ 个品牌中选择一些商品做营销活动. 已知每个商品的营销成本, 商品预期的收益和营销的预算. 在总营销费用不超过预算且每个品牌选中商品数量有限制的前提下, 如何选择商品使得预期的收益最大化?

我们先考虑一个直观的数学模型.

indices

$i$ -- 品牌
$j$ -- 商品

parameters

$a_{i, j} \in \{0, 1\}$ -- 商品 $j$ 是否属于品牌 $i$
$p_j$ -- 商品 $j$ 的预期收益
$c_j$ -- 商品 $j$ 的营销成本
$b_i$ -- 选中品牌 $i$ 的商品数量上限
$d$ -- 营销费用的总预算

decision variables

$x_j\in \{0,1\}$ -- 是否选择商品 $j$

模型1

$\begin{aligned} \max\ & \sum_j p_j x_j \\ \text{s.t. } & \sum_j c_j x_j \leq d \\ & \sum_j a_{i,j} x_j \leq b_i, \quad \forall i\\ & x_j \in \{0, 1\}. \end{aligned}$

当品牌和商品数量较大时, 例如1000品牌和10万商品, 这时参数 $a_{i,j}$ 的规模是1亿, 因此直接求解非常困难. 注意到每个品牌的商品是不一样的, 因而矩阵 $A = (a_{i,j})$ 非常稀疏, 我们可以把每个品牌中的商品分开考虑, 得到如下模型.

indices

$i$ -- 品牌
$j$ -- 商品

parameters

$n_i$ -- 品牌 $i$ 的商品数量
$p_{i,j}$ -- 品牌 $i$ 中商品 $j$ 的预期收益, $j=1, 2, \ldots, n_i$
$c_{i,j}$ -- 品牌 $i$ 中商品 $j$ 的营销成本, $j=1, 2, \ldots, n_i$
$b_i$ -- 品牌 $i$ 可以被选中的商品数量上限
$d$ -- 营销费用的总预算

decision variables

$x_{i,j}\in \{0, 1\}$ -- 是否选择品牌 $i$ 中的商品 $j$ , $j=1, 2, \ldots, n_i$

模型2
$\begin{aligned} \max\ & \sum_{i=1}^m\sum_{j=1}^{n_i} p_{i,j} x_{i,j} \\ \text{s.t. } & \sum_{i=1}^m\sum_{j=1}^{n_i} c_{i,j} x_{i,j} \leq d \\ & \sum_{j=1}^{n_i} x_{i,j} \leq b_i, \quad \forall i\\ & x_{i,j} \in \{0, 1\}. \end{aligned}$

模型1和模型2本质上是相同的, 因此当品牌数和商品数量非常大时直接求解模型2依然非常困难. 下面我们使用DW分解进行求解.

令 $P_i = \{x_{i,j} \mid \sum_{j=1}^{n_i} x_{i,j} \leq b_i \}$ , $i=1,2, \ldots, m$ . 注意到 $P_i$ 是有界的, 我们用
$v_{i,1}^{(k)}, v_{i,2}^{(k)}, \ldots, v_{i,p_i}^{(k)}, \quad k=1, 2, \ldots, p_i$

代表 $P_i$ 的顶点, 因此 $x_{i,j}$ 可以被表示成如下形式:
$x_{i,j} = \sum_{k=1}^{p_i}\lambda_{i,k}\cdot v_{i,j}^{(k)}, \quad \text{其中 } \sum_{k=1}^{p_i} \lambda_{i,k} = 1,\ \lambda_{i,k}\geq 0.$
把它代入模型2中我们得到主问题的形式.

主问题
$\begin{aligned} \max\ & \sum_{i=1}^m\sum_{j=1}^{n_i} \sum_{k=1}^{p_i}\lambda_{i,k} (p_{i,j}v_{i,j}^{(k)}) \\ \text{s.t. } & \sum_{i=1}^m \sum_{j=1}^{n_i} \sum_{k=1}^{p_i}\lambda_{i,k}(c_{i,j}v_{i,j}^{(k)}) \leq d \\ & \sum_{k=1}^{p_i} \lambda_{i,k} =1, \quad \forall i \\ & \lambda_{i,k} \geq 0, \quad \forall i, j, k. \end{aligned}$

定义对偶变量 $y$ 和 $z_i$ , $i=1,2,\ldots, m$ . 计算 $\lambda_{i,k}$ 对应的reduced cost:
$\alpha_i^{(k)} = \sum_{j=1}^{n_i} p_{i,j}v_{i,j}^{(k)} - \sum_{j=1}^{n_i}y\cdot c_{i,j}v_{i,j}^{(k)} - z_i.$

注意: 主问题是最大化问题, 因此 $\alpha_i^{(k)}>0$ 意味着可以提升主问题的目标函数值. 我们有:

子问题是最大化问题.
当所有为 $\alpha_i^{(k)} \leq 0$ 时主问题达到最优.

子问题 - $i$
$\begin{aligned} \min\ & \sum_{j=1}^{n_i}( p_{i,j}-yc_{i,j})x_{i,j} \\ \text{s.t. } & \sum_{j=1}^{n_i}x_{i,j} \leq b_i \\ & x_{i,j} \in \{0, 1\}. \end{aligned}$

初始化

对任意的 $i=1, 2, \ldots, m$ , 定义向量:
$v_i = (0, 0, \ldots, 0)^T \in \mathbb{R}^{b_i}.$
显然 $v_i$ 是每个约束 $i$ 的可行解, 即 $\sum_{j=1}^{b_i} v_{i,j} \leq b_i$ , $\forall i$ . 我们用 $v_1, v_2, \ldots, v_m$ 作为主问题初始化的顶点.

Remark. 前文的推导过程要求 $v_i$ 是多面体的顶点, 但上面 $v_i$ 的定义并不满足此条件. 这么做可行的原因是任意可行解本身就是多面体顶点的凸组合.

求解

求解的基本步骤如下:

初始化主问题, 求解子问题的输入参数 $y$
求解 $m$ 个子问题，分别计算 $\lambda_{i,k}$ 对应的Reduced Cost $\alpha_i^{(k)}$ . 如果 $\alpha_i^{(k)}>0$ , 则把对应的解 $v_i^{(k)}$ 加入到主问题. ( $k$ 可以理解为迭代的次数)
如果所有的 $\alpha_i^{(k)} < 0$ , 则停止迭代；否则迭代求解主问题和子问题直到满足停止条件.

Python实现

主问题模型

# master_model.py

from ortools.linear_solver import pywraplp


class MasterModel(object):

    def __init__(self, p, v, c, d):
        """
        :param p: p[i][j]代表品牌i中商品j的预期收益
        :param v: v[i]代表第i个子问题的解
        :param c: c[i][j]代表品牌i中商品j的营销成本
        :param d: scalar, 总预算
        """
        self._solver = pywraplp.Solver('MasterModel',
                                       pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
        self._p = p
        self._v = v
        self._c = c
        self._d = d
        self._la = None  # 决策变量lambda
        self._constraint_y = None  # 约束
        self._constraint_z = []  # 约束
        self._solution_la = None  # 计算结果

    def _init_decision_variables(self):
        self._la = [[]] * len(self._v)
        self._solution_la = [[]] * len(self._v)  # 初始化保存结果的变量
        for i in range(len(self._v)):
            self._la[i] = [[]] * len(self._v[i])
            self._solution_la[i] = [[]] * len(self._v[i])  # 初始化保存结果的变量
            for k in range(len(self._v[i])):
                self._la[i][k] = self._solver.NumVar(0, 1, 'la[%d][%d]' % (i, k))

    def _init_constraints(self):
        self._constraint_y = self._solver.Constraint(0, self._d)
        for i in range(len(self._v)):
            for k in range(len(self._v[i])):
                f = 0
                for j in range(len(self._v[i][k])):
                    f += self._c[i][j] * self._v[i][k][j]
                self._constraint_y.SetCoefficient(self._la[i][k], f)

        self._constraint_z = [None] * len(self._v)
        for i in range(len(self._v)):
            self._constraint_z[i] = self._solver.Constraint(1, 1)
            for k in range(len(self._la[i])):
                self._constraint_z[i].SetCoefficient(self._la[i][k], 1)

    def _init_objective(self):
        obj = self._solver.Objective()
        for i in range(len(self._v)):
            for k in range(len(self._v[i])):
                f = 0
                for j in range(len(self._v[i][k])):
                    f += self._p[i][j] * self._v[i][k][j]
                obj.SetCoefficient(self._la[i][k], f)
        obj.SetMaximization()

    def solve(self):
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()
        self._solver.Solve()
        # 保存计算结果
        for i in range(len(self._v)):
            for k in range(len(self._v[i])):
                self._solution_la[i][k] = self._la[i][k].solution_value()

    def get_solution_value(self):
        return self._solution_la

    def get_y(self):
        """ 获取对偶变量y的值
        """
        return self._constraint_y.dual_value()

    def get_zi(self, i):
        """ 获取对偶变量z[i]的值
        """
        return self._constraint_z[i].dual_value()

    def get_obj_value(self):
        res = 0
        for i in range(len(self._p)):
            for k in range(len(self._v[i])):
                for j in range(len(self._p[i])):
                    res += self._solution_la[i][k] * self._p[i][j] * self._v[i][k][j]
        return res

    def get_solution_x(self):
        """ 得到原问题的解.  x[i][j] = sum(la[i][k] * v[i][k][j]) over k.
        """

        x = [[]] * len(self._v)
        for i in range(len(self._v)):
            x[i] = [0] * len(self._v[i][0])

        for i in range(len(self._v)):
            for k in range(len(self._v[i])):
                for j in range(len(self._v[i][k])):
                    x[i][j] += self._solution_la[i][k] * self._v[i][k][j]
        return x

子问题模型

# sub_model.py

from ortools.linear_solver import pywraplp
import numpy as np


class SubModel(object):
    """ 子问题i
    """
    def __init__(self, pi, ci, y, bi):
        """ 下标i忽略
        :param pi: list, pi := p[i] = [p1, p2, ..., ]
        :param ci: list, ci := c[i] = [c1, c2, ..., ]
        :param y: scalar
        :param bi: scalar
        """
        self._solver = pywraplp.Solver('MasterModel',
                                       pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
        self._pi = pi
        self._ci = ci
        self._y = y
        self._bi = bi
        self._x = None  # 决策变量
        self._solution_x = None  # 计算结果

    def _init_decision_variables(self):
        self._x = [None] * len(self._pi)
        for j in range(len(self._pi)):
            self._x[j] = self._solver.IntVar(0, 1, 'x[%d]' % j)

    def _init_constraints(self):
        ct = self._solver.Constraint(0, self._bi)
        for j in range(len(self._pi)):
            ct.SetCoefficient(self._x[j], 1)

    def _init_objective(self):
        obj = self._solver.Objective()
        for j in range(len(self._pi)):
            obj.SetCoefficient(self._x[j], self._pi[j] - self._y * self._ci[j])
        obj.SetMaximization()

    def solve(self):
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()
        self._solver.Solve()
        self._solution_x = [s.solution_value() for s in self._x]

    def get_solution_x(self):
        return self._solution_x

    def get_obj_value(self):
        p = np.array(self._pi)
        c = np.array(self._ci)
        x = np.array(self._solution_x)
        return sum((p - self._y * c) * x)

DW分解的求解过程

# dw_proc.py

from master_model import MasterModel
from sub_model import SubModel


class DWProc(object):

    def __init__(self, p, c, d, b, max_iter=1000):
        """
        :param p: p[i][j]代表品牌i中商品j的预期收益
        :param c: c[i][j]代表品牌i中商品j的营销成本
        :param d: 总营销成本, int
        :param b: b[i]代表选中品牌i的商品数量限制
        """
        self._p = p
        self._c = c
        self._d = d
        self._b = b
        self._v = None  # 待初始化
        self._max_iter = max_iter
        self._iter_times = 0
        self._status = -1
        self._reduced_costs = [1] * len(self._p)
        self._solution_x = None  # 计算结果
        self._obj_value = 0  # 目标函数值

    def _stop_criteria_is_satisfied(self):
        """ 根据reduced cost判断是否应该停止迭代.
        注意: 这是最大化问题, 因此所有子问题对应的reduced cost <= 0时停止.
        """
        st = [0] * len(self._reduced_costs)
        for i in range(len(self._reduced_costs)):
            if self._reduced_costs[i] < 1e-6:
                st[i] = 1
        if sum(st) == len(st):
            self._status = 0
            return True
        if self._iter_times >= self._max_iter:
            if self._status == -1:
                self._status = 1
            return True
        return False

    def _init_v(self):
        """ 初始化主问题的输入
        """
        self._v = [[]] * len(self._p)
        for i in range(len(self._p)):
            self._v[i] = [[0] * len(self._p[i])]

    def _append_v(self, i, x):
        """ 把子问题i的解加入到主问题中

        :param x: 子问题i的解
        """
        self._v[i].append(x)

    def solve(self):
        # 初始化主问题并求解
        self._init_v()
        mp = MasterModel(self._p, self._v, self._c, self._d)
        mp.solve()
        self._iter_times += 1
        # 迭代求解主问题和子问题直到满足停止条件
        while not self._stop_criteria_is_satisfied():
            # 求解子问题
            print("==== iter %d ====" % self._iter_times)
            for i in range(len(self._p)):
                # 求解子问题
                sm = SubModel(self._p[i], self._c[i], mp.get_y(), self._b[i])
                sm.solve()
                # 更新reduced cost
                self._reduced_costs[i] = sm.get_obj_value() - mp.get_zi(i)
                # 把子问题中满足条件的解加入到主问题中
                if self._reduced_costs[i] > 0:
                    self._append_v(i, sm.get_solution_x())
                print(">>> Solve sub problem %d, reduced cost = %f" % (i, self._reduced_costs[i]))

            # 求解主问题
            mp = MasterModel(self._p, self._v, self._c, self._d)
            mp.solve()

            self._iter_times += 1

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

    def print_info(self):
        print("==== Result Info  ====")
        print(">>> objective value =", self._obj_value)
        print(">>> Solution")
        sku_list = [[]] * len(self._solution_x)
        for i in range(len(self._solution_x)):
            sku_list[i] = [j for j in range(len(self._solution_x[i])) if self._solution_x[i][j] > 0]
        for i in range(len(self._solution_x)):
            print("brand %d, sku list:" % i, sku_list[i])

主函数

# main.py

from data import p, c, b, d  # instance data
from dw_proc import DWProc


if __name__ == '__main__':
    dw = DWProc(p, c, d, b)
    dw.solve()
    dw.print_info()

完整代码

参考文献

George B. Dantzig; Philip Wolfe. Decomposition Principle for Linear Programs. Operations Research. Vol 8: 101–111, 1960. ↩

线性规划技巧: Dantzig&Wolfe Decomposti

问题描述

主问题

子问题

例子: 一个选品问题

Python实现

参考文献

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读