线性规划技巧: Dantzig&Wolfe Decompostion

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

问题描述

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

$$\begin{array}{rlrlrlrlr}\min & c_1^T{x}_1& +& c_2^T{x}_2& +& \cdots & +& c_m^T{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\ge 0,& \forall j& & & & & & \end{array}$$

其中$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}$.

说明

1. 上述规划共$m+1$组约束. 第一组约束包含了所有的决策变量, 称为连接约束(Linking constraints), 接下来是$m$组独立的约束.
2. $c_j,x_j,b_i$都是向量.
3. 上述规划如果写成标准形式 min ⁡ x { c T x ∣ A x = b , x ≥ 0 } \min_x \{c^Tx | Ax = b, x\geq 0\} minx​{cTx∣Ax=b,x≥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 ∈ R n i ∣ A i , i x i = b i } P_i = \{ x_i \in \mathbb{R}^{n_i} \mid A_{i,i}x_i = b_i\} Pi​={xi​∈Rni​∣Ai,i​xi​=bi​}. 如果$P_i$是有界的, 那么$x_i$可以表示成$P_i$顶点的凸组合. 设$v_{i,1},v_{i,2},\dots ,v_{i,p_i}$代表$P_i$的顶点(Extreme point), 则存在${\lambda }_{i,j}\ge 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},\dots ,w_{i,q_i}$代表极射线(Extreme ray). 则存在${\lambda }_{i,j}\ge 0$且${\sum }_{j=1}^{p_i}{\lambda }_{i,j}=1$, ${\mu }_{i,j}\ge 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,\dots ,P_m$所有的顶点和极射线. 把上式代入原问题得到主问题(Master problem):

主问题

$$\begin{array}{rl}\min & \sum _{i=1}^m\sum _{j=1}^{p_i}{\lambda }_{i,j}(c_i^T{v}_{i,j})+\sum _{i=1}^m\sum _{j=1}^{q_i}{\mu }_{i,j}(c_i^T{w}_{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,\forall i\\ & {\lambda }_{i,j},{\mu }_{i,j}\ge 0,\forall i,j.\end{array}$$

说明

1. ${\lambda }_{i,j}$, ${\mu }_{i,j}$是决策变量.
2. 约束的数量为$m_0+m$(第一个等式包含了$m_0$个约束). 原问题的约束数量是${\sum }_{i=0}^m{m}_i$. 因此主问题的约束数量明显减少了.
3. 但是主问题的变量显著增加了(${\lambda }_{i,j},{\mu }_{i,j}$对应所有的顶点和极射线). 与列生成技巧相似, 主问题一开始只考虑两个可行解(对应顶点和极射线). 通过求解子问题得到新的顶点或极射线.
4. 在实际问题中, 一般情况下$P_i$都是有界的. 此时主问题可以简化成如下形式.
   $$\begin{array}{rl}\min & \sum _{i=1}^m\sum _{j=1}^{p_i}{\lambda }_{i,j}(c_i^T{v}_{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,\forall i\\ & {\lambda }_{i,j}\ge 0,\forall i,j.\end{array}$$

子问题

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

$$\begin{array}{rl}& {\alpha }_{i,j}=c_i^T{v}_{i,j}-y^T{A}_{0,i}{v}_{i,j}-z_i\\ & {\beta }_{i,j}=c_i^T{w}_{i,j}-y^T{A}_{0,i}{w}_{i,j}\end{array}$$

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

子问题 - $i$

$$\begin{array}{rl}\min & f_i:=(c_i-A_{0,i}^{T}y{)}^T{x}_i\\ \text{s.t. }& A_{i,i}{x}_i=b_i\\ & x_i\ge 0.\end{array}$$

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

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

例子: 一个选品问题

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

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

indices

• $i$ – 品牌
• $j$ – 商品

parameters

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

decision variables

• x j ∈ { 0 , 1 } x_j\in \{0,1\} xj​∈{0,1} – 是否选择商品$j$

模型1

max ⁡   ∑ j p j x j s.t.  ∑ j c j x j ≤ d ∑ j a i , j x j ≤ b i , ∀ i x j ∈ { 0 , 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} max s.t. ​j∑​pj​xj​j∑​cj​xj​≤dj∑​ai,j​xj​≤bi​,∀ixj​∈{0,1}.​

当品牌和商品数量较大时, 例如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,\dots ,n_i$
• $c_{i,j}$ – 品牌$i$中商品$j$的营销成本, $j=1,2,\dots ,n_i$
• $b_i$ – 品牌$i$可以被选中的商品数量上限
• $d$ – 营销费用的总预算

decision variables

• x i , j ∈ { 0 , 1 } x_{i,j}\in \{0, 1\} xi,j​∈{0,1} – 是否选择品牌$i$中的商品$j$, $j=1,2,\dots ,n_i$

模型2
max ⁡   ∑ i = 1 m ∑ j = 1 n i p i , j x i , j s.t.  ∑ i = 1 m ∑ j = 1 n i c i , j x i , j ≤ d ∑ j = 1 n i x i , j ≤ b i , ∀ i x i , j ∈ { 0 , 1 } . \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} max s.t. ​i=1∑m​j=1∑ni​​pi,j​xi,j​i=1∑m​j=1∑ni​​ci,j​xi,j​≤dj=1∑ni​​xi,j​≤bi​,∀ixi,j​∈{0,1}.​

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

令 P i = { x i , j ∣ ∑ j = 1 n i x i , j ≤ b i } P_i = \{x_{i,j} \mid \sum_{j=1}^{n_i} x_{i,j} \leq b_i \} Pi​={xi,j​∣∑j=1ni​​xi,j​≤bi​}, $i=1,2,\dots ,m$. 注意到$P_i$是有界的, 我们用
$$v_{i,1}^{(k)},v_{i,2}^{(k)},\dots ,v_{i,p_i}^{(k)},k=1,2,\dots ,p_i$$

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

主问题
$$\begin{array}{rl}\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)})\le d\\ & \sum _{k=1}^{p_i}{\lambda }_{i,k}=1,\forall i\\ & {\lambda }_{i,k}\ge 0,\forall i,j,k.\end{array}$$

定义对偶变量$y$和$z_i$, $i=1,2,\dots ,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$意味着可以提升主问题的目标函数值. 我们有:

1. 子问题是最大化问题.
2. 当所有为${\alpha }_i^{(k)}\le 0$时主问题达到最优.

子问题 - $i$
min ⁡   ∑ j = 1 n i ( p i , j − y c i , j ) x i , j s.t.  ∑ j = 1 n i x i , j ≤ b i x i , j ∈ { 0 , 1 } . \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} min s.t. ​j=1∑ni​​(pi,j​−yci,j​)xi,j​j=1∑ni​​xi,j​≤bi​xi,j​∈{0,1}.​

初始化

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

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

求解

求解的基本步骤如下:

1. 初始化主问题, 求解子问题的输入参数$y$
2. 求解$m$个子问题，分别计算${\lambda }_{i,k}$对应的Reduced Cost ${\alpha }_i^{(k)}$. 如果${\alpha }_i^{(k)}>0$, 则把对应的解$v_i^{(k)}$加入到主问题. ($k$可以理解为迭代的次数)
3. 如果所有的${\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()

完整代码

参考文献

---

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