旅行商問題（TSP）系列(1)

嘗試使用各種優(yōu)化算法诗箍，對(duì)旅行商問題（TSP）及多旅行商問題（MTSP）進(jìn)行解決卷中。主要使用語言為python。

題目

(1) 在地圖上給定30個(gè)節(jié)點(diǎn)的位置褥符，標(biāo)號(hào)為0--29，每兩個(gè)節(jié)點(diǎn)之間都有一條路徑相互連通抚垃，求一條回路能遍歷每一個(gè)節(jié)點(diǎn)，并且令其總長度盡可能小趟大。

(2) 在第一問的基礎(chǔ)上鹤树，求四條回路，它們加起來能遍歷每一個(gè)節(jié)點(diǎn)逊朽，并且令這四條回路的總長度盡可能小罕伯。

(3) 在第二問的基礎(chǔ)上，增加限制條件叽讳，要求這四條回路都經(jīng)過0節(jié)點(diǎn)追他，令這四條回路的總長度盡可能小。

附件：40個(gè)節(jié)點(diǎn)的位置（按順序分別為0節(jié)點(diǎn)岛蚤，1節(jié)點(diǎn)邑狸，……29節(jié)點(diǎn)）：

(45, 0), (67, 17), (22, 56), (42, 9), (48, 13),
(4, 66), (58, 9), (5, 0), (81, 41), (97, 21),
(63, 61), (27, 76), (46, 73), (41, 63), (19, 83),
(86, 46), (73, 83), (9, 25), (88, 96), (91, 63),
(4, 76), (51, 56), (62, 71), (76, 63), (90, 40),
(86, 15), (25, 63), (15, 54), (18, 37), (61, 10)

30個(gè)點(diǎn)的位置 可視化

數(shù)據(jù)預(yù)處理

先對(duì)點(diǎn)的位置數(shù)據(jù)進(jìn)行預(yù)處理，并準(zhǔn)備幾個(gè)常用函數(shù)涤妒。

# basic.py
# 對(duì)點(diǎn)的位置數(shù)據(jù)進(jìn)行預(yù)處理
# 以及后面會(huì)用到的常用函數(shù)
from math import sqrt
import numpy as np
import matplotlib.pyplot as plt

# 節(jié)點(diǎn)數(shù)量
num_pos = 30

# 各節(jié)點(diǎn)位置
pos = (
    (45, 0), (67, 17), (22, 56), (42, 9), (48, 13), 
    (4, 66), (58, 9), (5, 0), (81, 41), (97, 21), 
    (63, 61), (27, 76), (46, 73), (41, 63), (19, 83), 
    (86, 46), (73, 83), (9, 25), (88, 96), (91, 63), 
    (4, 76), (51, 56), (62, 71), (76, 63), (90, 40), 
    (86, 15), (25, 63), (15, 54), (18, 37), (61, 10)
    )

# 計(jì)算兩點(diǎn)間距離
def get_dist(pos_1, pos_2):
    return sqrt((pos_1[0] - pos_2[0])**2 + (pos_1[1] - pos_2[1])**2)

# 各節(jié)點(diǎn)之間的距離
dist = [[get_dist(pos_1, pos_2) for pos_2 in pos] for pos_1 in pos]


# 繪制散點(diǎn)圖
def show_points(path):
    # 例：{v_0, v_2, v_1, v_3}: path = [0, 2, 1, 3]
    x = [pos[v_i][0] for v_i in path]
    y = [pos[v_i][1] for v_i in path]
    # 繪制散點(diǎn)圖
    plt.scatter(x, y)
    # 標(biāo)號(hào)
    for i in range(len(path)):
        plt.annotate(path[i], (x[i] + 1, y[i] - 1))
    # plt.show()


# 繪制回路圖
def show_path(path):
    # 例：v_0-v_2-v_1-v_3-v_0: path = [0, 2, 1, 3, 0]
    x = [pos[v_i][0] for v_i in path]
    y = [pos[v_i][1] for v_i in path]
    # 繪制散點(diǎn)圖
    plt.scatter(x, y)
    # 繪制連線
    plt.plot(x, y)
    # 標(biāo)號(hào)
    for i in range(len(path)):
        plt.annotate(path[i], (x[i] + 1, y[i] - 1))
    # plt.show()

第一問：狀壓dp

思路

(1) 在地圖上給定30個(gè)節(jié)點(diǎn)的位置单雾，標(biāo)號(hào)為0--29，每兩個(gè)節(jié)點(diǎn)之間都有一條路徑相互連通她紫，求一條回路能遍歷每一個(gè)節(jié)點(diǎn)硅堆，并且令其總長度盡可能小。

很顯然贿讹，第一問便是求最短哈密頓回路渐逃，也就是很經(jīng)典的旅行商問題。旅行商問題的可行解是所有頂點(diǎn)的全排列民褂，隨著頂點(diǎn)數(shù)的增加茄菊，會(huì)產(chǎn)生組合爆炸，難以用窮舉法解出答案赊堪。[1]但考慮到這道題只有30個(gè)點(diǎn)买羞，并不算太多（后面被計(jì)算量狠狠打臉），而且也不確定是否能接受非最優(yōu)解雹食，我一開始決定先用搜索法窮舉所有可能畜普，找出全局最優(yōu)解。

在方法選擇上群叶，深度優(yōu)先搜索+剪枝顯然不夠給力吃挑，于是我選擇了狀壓dp（狀態(tài)壓縮動(dòng)態(tài)規(guī)劃）钝荡。[2]

用離散數(shù)學(xué)的知識(shí)，對(duì)這道題進(jìn)行抽象化描述：

已知無向完全圖G=<V,E>中舶衬，V=\{v_0, v_1, v_2, ..., v_{29}\}埠通。設(shè)狀態(tài)s為一個(gè)長度為30的01串，用來表示V的一個(gè)子集V_s逛犹。比如用s_{11}=000000000000000000000000001101端辱，表示V的一個(gè)子集V_{11} = \{v_0, v_2, v_3\}。通過這種二進(jìn)制的表示方法虽画，我們可以把30個(gè)節(jié)點(diǎn)的所有狀態(tài)（子集）壓縮到一個(gè)int里舞蔽，此時(shí)s的取值范圍是[0, 2^{30} -1]。

設(shè)dist(i, j)為v_i和v_j間的距離码撰。設(shè)min\_length(s, i)為：G的完全無向子圖G_s=<V_s,E_s>中渗柿，以v_0為起點(diǎn)，v_i為終點(diǎn)的哈密頓通路（經(jīng)過s中所有節(jié)點(diǎn)一次且僅一次的通路）的最短長度（v_0,v_i \in s）脖岛。

若min\_length(s +v_i, i) \ (v_i \notin s)所指的這條哈密頓通路為v_0 ... v_jv_i朵栖，則d(v_0...v_jv_i) = d(v_0...v_j) + d(v_iv_j)，因而我們可以得到狀態(tài)轉(zhuǎn)移方程：
min\_length(s + v_i, i) = \min_{j \in s, \ j \neq 0}{\left(min\_length(s, j) + dist(i, j) \right)} \\ (v_0 \in s, v_i \notin s)
該狀態(tài)轉(zhuǎn)移方程的含義是：在到達(dá)v_i之前柴梆，我們先到達(dá)的是v_j陨溅；在到達(dá)s+v_i這個(gè)狀態(tài)前，我們先到達(dá)的是s這個(gè)狀態(tài)绍在。因此只要得到所有min\_length(s, j)的值声登，我們就能得知min\_length(s + v_i, i)的值。換言之揣苏，每次到達(dá)一個(gè)min\_length(s, j)時(shí)悯嗓，我們都可以更新所有min\_length(s + v_i, i)的值。(j \in s, i \notin s)

我們發(fā)現(xiàn)卸察，壓縮到int中時(shí)脯厨，s + v_i記作s | (1 << i)，換言之int(s + v_i) > int(s)恒成立坑质。這啟示我們可以以一個(gè)有序（從小到大）的順序遍歷s合武，使得每次抵達(dá)一個(gè)s+v_i時(shí)，都可以保證它已經(jīng)被所有的s更新過涡扼。這樣稼跳，我們就可以通過動(dòng)態(tài)規(guī)劃，得到最終G中以v_0為起點(diǎn)吃沪，每一個(gè)以v_i為終點(diǎn)的哈密頓通路的最短長度汤善。那么，可以通過d(v_0..v_i v_0) = d(v_0...v_i) + d(v_i v_0)，得到最短哈密頓回路的長度ans\_length：
ans\_length = \min_{i \in s, \ i \neq 0}{\left( min\_length(s, i) + dist(i, 0) \right)}

代碼

這種算法用c++跑肯定會(huì)快一些红淡，但是python比較方便（可以作圖）不狮。

# q_1_1.py
# 狀壓dp求解最短哈密頓回路
from basic import *

# 先用較少的節(jié)點(diǎn)數(shù)做測試
num_pos = 10
# 無窮大常量
inf = 987654321
# s中以0為起點(diǎn)，i為終點(diǎn)的最短哈密頓通路長度
min_dist = [[inf for i in range(num_pos)] for s in range(1 << num_pos)]
# s中以0為起點(diǎn)在旱，i為終點(diǎn)的最短哈密頓通路
# 例: v_0-v_3-v_1-v_2表示為[3, 1, 2]
min_path = [[[] for i in range(num_pos)] for s in range(1 << num_pos)]

print('init succeed')


# 狀壓dp搜索
for s in range(1, 1 << num_pos, 2):
    # 考慮節(jié)點(diǎn)集合s必須包括節(jié)點(diǎn)0
    if not (s & 1):
        continue
    for j in range(1, num_pos):
        # 終點(diǎn)i需在當(dāng)前考慮節(jié)點(diǎn)集合s內(nèi)
        if not (s & (1 << j)):
            continue
        if s == int((1 << j) | 1):
            # 若考慮節(jié)點(diǎn)集合s僅含節(jié)點(diǎn)0和節(jié)點(diǎn)j摇零，dp邊界，賦予初值
            # print('j:', j)
            min_path[s][j] = [j]
            min_dist[s][j] = dist[0][j]

        # 枚舉下一個(gè)節(jié)點(diǎn)i桶蝎，更新
        for i in range(1, num_pos):
            # 下一個(gè)節(jié)點(diǎn)i需在考慮節(jié)點(diǎn)集合s外
            if s & (1 << i):
                continue
            # 更新min_dist[s + i][i], min_path[s + i][i]
            if min_dist[s][j] + dist[j][i] < min_dist[s | (1 << i)][i]:
                min_path[s | (1 << i)][i] = min_path[s][j] + [i]
                min_dist[s | (1 << i)][i] = min_dist[s][j] + dist[j][i]


ans_dist = inf
ans_path = [] 
# 求最終最短哈密頓回路
for i in range(1, num_pos):
    if min_dist[(1 << num_pos) - 1][i] + dist[i][0] < ans_dist:
        # 更新驻仅，回路化
        ans_path = [0] + min_path[s][i] + [0]
        ans_dist = min_dist[(1 << num_pos) - 1][i] + dist[i][0]

# 輸出結(jié)果
print()
print('min length:', ans_dist)
print('path:', ans_path)

show_path(ans_path)
plt.show()

運(yùn)行結(jié)果

在第5行注意到，我只使用10個(gè)節(jié)點(diǎn)做測試登渣，結(jié)果如下：

1595563051243.png
Figure_2-1595560139309.png

可以看到10節(jié)點(diǎn)時(shí)噪服，狀壓dp跑的很快，解的形狀看上去也完全沒問題绍豁。接下來慢慢放開節(jié)點(diǎn)數(shù)量，看看跑得如何牙捉。

節(jié)點(diǎn)數(shù) 所用時(shí)間 最短長度
10 0.8s 282.9
15 2.1s 325.8
18 14.4s 366.0
20 100.3s 401.3
21 144.0s 405.1
22 1013.3s 409.7
... ... ...

……我佛了竹揍。

（這個(gè)結(jié)論或許可以記一下？22個(gè)點(diǎn)以內(nèi)的TSP都可以用狀壓dp找最優(yōu)解）

非常明顯邪铲，這個(gè)算法的時(shí)間開銷是我們不能接受的（想打死之前覺得30不算多的我）芬位。掐指一算，該算法的時(shí)間復(fù)雜度是O(n^2 \cdot 2^n)带到，空間復(fù)雜度是O(n \cdot 2^n)昧碉，（光是數(shù)組初始化也要很長時(shí)間）。這個(gè)令人絕望的增長速度也打消了我“用c++重跑一遍”的念頭揽惹。寫到這里被饿，我感覺30個(gè)點(diǎn)的最優(yōu)解求不出來，決定放棄狀壓dp搪搏，改用一些比較隨機(jī)化的算法求一個(gè)較優(yōu)解狭握。

第一問：模擬退火算法[3]

簡介

模擬退火算法是一種啟發(fā)式搜索算法，即按照預(yù)定的控制策略進(jìn)行搜索疯溺，在搜索過程中獲取的中間信息將用來改進(jìn)控制策略 论颅。它的思想借鑒于固體的退火原理，當(dāng)固體的溫度很高的時(shí)候囱嫩，內(nèi)能比較大恃疯，固體的內(nèi)部粒子處于快速無序運(yùn)動(dòng)，當(dāng)溫度慢慢降低的過程中墨闲，固體的內(nèi)能減小今妄，粒子的慢慢趨于有序，最終，當(dāng)固體處于常溫時(shí)蛙奖，內(nèi)能達(dá)到最小潘酗，此時(shí)，粒子最為穩(wěn)定雁仲。模擬退火算法便是基于這樣的原理設(shè)計(jì)而成仔夺。

模擬退火算法從某一高溫出發(fā)，在高溫狀態(tài)下計(jì)算初始解攒砖，然后以預(yù)設(shè)的鄰域函數(shù)產(chǎn)生一個(gè)擾動(dòng)量缸兔，從而得到新的狀態(tài)，即模擬粒子的無序運(yùn)動(dòng)吹艇，比較新舊狀態(tài)下的能量惰蜜，即目標(biāo)函數(shù)的解。如果新狀態(tài)的能量小于舊狀態(tài)受神，則狀態(tài)發(fā)生轉(zhuǎn)化抛猖；如果新狀態(tài)的能量大于舊狀態(tài)，則以一定的概率準(zhǔn)則發(fā)生轉(zhuǎn)化鼻听。當(dāng)狀態(tài)穩(wěn)定后财著，便可以看作達(dá)到了當(dāng)前狀態(tài)的最優(yōu)解，便可以開始降溫撑碴，在下一個(gè)溫度繼續(xù)迭代撑教，最終達(dá)到低溫的穩(wěn)定狀態(tài)，便得到了模擬退火算法產(chǎn)生的結(jié)果醉拓。

在這道題中伟姐，將回路\Gamma = v_{i_0} v_{i_1} v_{i_2} ... v_{i_{29}} v_{i_0}表示為x = [i_0, i_1, i_2, ..., i_29]，將x在附近鄰域隨機(jī)擾動(dòng)設(shè)為令x中多個(gè)元素隨機(jī)交換位置亿卤，便可以使x的取值范圍覆蓋整個(gè)解空間愤兵。

代碼

于是我便得到一個(gè)（比較粗糙的）模擬退火算法的代碼。

# q_1_1.py
# 模擬退火算法求解最短哈密頓回路
from basic import *
from random import randint, uniform, shuffle
from math import exp, log

# 計(jì)算回路長度
def get_path_length(x):
    length = 0
    for i in range(len(x) - 1):
        length += dist[x[i]][x[i + 1]]
    length += dist[len(x) - 1][0]
    return length

# 給予x一個(gè)隨機(jī)位移
def get_next_x(x):
    x_next = x.copy()
    # 多次交換
    for i in range(int(T / T_0 * 8 + 2)):
        i_1, i_2 = randint(0, num_pos - 1), randint(0, num_pos - 1)
        x_next[i_1], x_next[i_2] = x_next[i_2], x_next[i_1]
    return x_next

# 采納概率
# 自己瞎搗鼓出來的
def get_accpectable_probability(y, y_next, T):
    return exp(-T_0 / T * (y_next - y) / y)


# # 溫度初值
T_0 = 100000
T = T_0
# 溫度終止值
T_min = 30
# x: 回路
# 例: v_0-v_1-v_2-v_0: x = [0, 1, 2]
# x初值隨機(jī)
x = list(range(num_pos))
shuffle(x)
# y: 回路長度
y = get_path_length(x)
# 內(nèi)循環(huán)次數(shù)
k = 100
# 計(jì)數(shù)器
t = 0           

# 存儲(chǔ)求解過程
x_list = [x]
y_list = [y]

x_best = x.copy()
y_best = get_path_length(x)

# 開始模擬退火
while T > T_min:
    # 內(nèi)循環(huán)
    for i in range(k):
        y = get_path_length(x)

        # 給予x隨機(jī)位移
        x_next = get_next_x(x)

        # 試探新解
        y_next = get_path_length(x_next)
        if y_next < y:
            # 更優(yōu)解排吴，一定接受新解
            x = x_next
            # 更新最優(yōu)記錄
            if y_next < y_best:
                y_best = y_next
                x_best = x_next.copy()
        else:
            # 更劣解恐似，一定概率接受新解
            p = get_accpectable_probability(y, y_next, T)
            r = uniform(0, 1)
            if r < p:
                x = x_next

    # 做記錄
    x_list.append(x)
    y_list.append(y_next)

    t += 1
    if t % 1000 == 0:
        print(T)

    # 降溫
    T = T_0 / (1 + t) # 快速降溫
    # T = T_0 / log(1 + t) # 經(jīng)典降溫，慢的一批傍念，令人暴躁


x_best += [x_best[0]]
# 輸出解
print('min length:', y_best)
print('path:', x_best)

plt.figure(1)
ax1 = plt.subplot(1, 2, 1)
ax2 = plt.subplot(1, 2, 2)

# 繪制連線圖
plt.sca(ax1)
show_path(x_best)

# 繪制退火圖
plt.sca(ax2)
x = np.linspace(1, len(x_list), len(x_list))
y = np.linspace(1, len(x_list), len(x_list))
for i in range(len(x_list)):
    y[i] = y_list[i]
plt.plot(x, y)

plt.title('模擬退火進(jìn)程', fontproperties = 'SimHei')
plt.xlabel('遍歷次數(shù)', fontproperties = 'SimHei')
plt.ylabel('回路長度', fontproperties = 'SimHei')

plt.show()

運(yùn)行結(jié)果

哪怕是在30個(gè)節(jié)點(diǎn)的情況下矫夷，速度依然飛快。多次運(yùn)行憋槐，結(jié)果如下：

運(yùn)行次數(shù) 回路長度 運(yùn)行時(shí)間
1 514.5 9.9s
2 486.7 11.3s
3 486.8 10.3s
4 497.4 11.5s
5 463.0 10.5s
6 529.2 11.2s
7 458.5 10.4s
8 455.4 10.8s
... ... ...

一些搜索結(jié)果如下双藕，可以看出，有時(shí)候找到的解明顯不太行阳仔，有時(shí)候找到的解很舒服忧陪。

Figure_3_1.png
Figure_3_2.png
1595566034252.png

總的來說扣泊，從右圖的模擬退火進(jìn)程可以看出，模擬退火受初值影響較大嘶摊，最終的結(jié)果主要取決于一開始那一小段能降到多低延蟹，在后續(xù)大部分時(shí)間難以進(jìn)一步優(yōu)化。因此叶堆，每一次跑的結(jié)果差異顯著阱飘，有時(shí)很好，有時(shí)很爛虱颗，全看運(yùn)氣沥匈，不能做到“最終都收斂于相近的答案”。

跑了四五十遍忘渔，得到的最短回路長度為453.2高帖，結(jié)果如下：

Figure_3_4.png

憑感覺，這個(gè)路線已經(jīng)找不出什么明顯的優(yōu)化空間了畦粮，將這個(gè)路線作為第一問的解答足以令人滿意散址。

待續(xù)

旅行商問題（TSP）系列(2)

參考資料

  1. 旅行商問題 百度百科 https://baike.baidu.com/item/%E6%97%85%E8%A1%8C%E5%95%86%E9%97%AE%E9%A2%98/7737042 ?

  2. 最短哈密頓路徑（位運(yùn)算+dp） solego https://blog.csdn.net/weixin_43900869/article/details/102934460 ?

  3. 模擬退火算法與其python實(shí)現(xiàn) WFRainn https://blog.csdn.net/wfrainn/article/details/80303138 ?

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末，一起剝皮案震驚了整個(gè)濱河市宣赔，隨后出現(xiàn)的幾起案子预麸，更是在濱河造成了極大的恐慌，老刑警劉巖拉背，帶你破解...
    沈念sama閱讀 222,590評(píng)論 6贊 517
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件师崎，死亡現(xiàn)場離奇詭異默终，居然都是意外死亡椅棺，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,157評(píng)論 3贊 399
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門齐蔽，熙熙樓的掌柜王于貴愁眉苦臉地迎上來两疚，“玉大人，你說我怎么就攤上這事含滴∮詹常” “怎么了？”我有些...
    開封第一講書人閱讀 169,301評(píng)論 0贊 362
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵谈况，是天一觀的道長勺美。 經(jīng)常有香客問我，道長碑韵，這世上最難降的妖魔是什么赡茸？ 我笑而不...
    開封第一講書人閱讀 60,078評(píng)論 1贊 300
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮祝闻，結(jié)果婚禮上占卧，老公的妹妹穿的比我還像新娘。我一直安慰自己，他們只是感情好华蜒，可當(dāng)我...
    茶點(diǎn)故事閱讀 69,082評(píng)論 6贊 398
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布辙纬。 她就那樣靜靜地躺著，像睡著了一般叭喜。 火紅的嫁衣襯著肌膚如雪贺拣。 梳的紋絲不亂的頭發(fā)上，一...
    開封第一講書人閱讀 52,682評(píng)論 1贊 312
  • 城市分裂傳說
    那天域滥，我揣著相機(jī)與錄音纵柿，去河邊找鬼。 笑死启绰，一個(gè)胖子當(dāng)著我的面吹牛昂儒，可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播委可，決...
    沈念sama閱讀 41,155評(píng)論 3贊 422
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼渊跋，長吁一口氣：“原來是場噩夢啊……” “哼！你這毒婦竟也來了着倾？” 一聲冷哼從身側(cè)響起拾酝，我...
    開封第一講書人閱讀 40,098評(píng)論 0贊 277
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對(duì)情侶失蹤，失蹤者是張志新（化名）和其女友劉穎卡者，沒想到半個(gè)月后蒿囤，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
    沈念sama閱讀 46,638評(píng)論 1贊 319
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡崇决，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 38,701評(píng)論 3贊 342
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年材诽，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片恒傻。...
    茶點(diǎn)故事閱讀 40,852評(píng)論 1贊 353
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡脸侥，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出盈厘，到底是詐尸還是另有隱情睁枕，我是刑警寧澤，帶...
    沈念sama閱讀 36,520評(píng)論 5贊 351
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布沸手，位于F島的核電站外遇，受9級(jí)特大地震影響，放射性物質(zhì)發(fā)生泄漏契吉。R本人自食惡果不足惜跳仿，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 42,181評(píng)論 3贊 335
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望栅隐。 院中可真熱鬧塔嬉，春花似錦玩徊、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,674評(píng)論 0贊 25
  • 一樁弒父案，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽读第。三九已至，卻和暖如春畔塔，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背鸯屿。 一陣腳步聲響...
    開封第一講書人閱讀 33,788評(píng)論 1贊 274
  • 情欲美人皮
    我被黑心中介騙來泰國打工澈吨， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留，地道東北人寄摆。 一個(gè)月前我還...
    沈念sama閱讀 49,279評(píng)論 3贊 379
  • 代替公主和親
    正文 我出身青樓谅辣，卻偏偏與公主長得像，于是被迫代替她去往敵國和親婶恼。 傳聞我的和親對(duì)象是個(gè)殘疾皇子桑阶，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,851評(píng)論 2贊 361