线性规划技巧: Dantzig&Wolfe Decompostion

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}.

说明

  1. 上述规划共m+1组约束. 第一组约束包含了所有的决策变量, 称为连接约束(Linking constraints), 接下来是m组独立的约束.
  2. c_j, x_j, b_i都是向量.
  3. 上述规划如果写成标准形式\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}

说明

  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{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}

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

  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} \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意味着可以提升主问题的目标函数值. 我们有:

  1. 子问题是最大化问题.
  2. 当所有为\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的定义并不满足此条件. 这么做可行的原因是任意可行解本身就是多面体顶点的凸组合.

求解

求解的基本步骤如下:

  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. ↩

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 205,236评论 6赞 478
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 87,867评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 151,715评论 0赞 340
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,899评论 1赞 278
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,895评论 5赞 368
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,733评论 1赞 283
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,085评论 3赞 399
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,722评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 43,025评论 1赞 300
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,696评论 2赞 323
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,816评论 1赞 333
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,447评论 4赞 322
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,057评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,009评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,254评论 1赞 260
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,204评论 2赞 352
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,561评论 2赞 343

推荐阅读更多精彩内容