999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于混合蟻群算法的重調度低碳城市配送

2022-10-17 13:53:36張明偉屈曉龍
計算機工程與設計 2022年10期
關鍵詞:優化

張明偉,李 波,屈曉龍

(1.天津仁愛學院 經濟與管理學院,天津 301636;2.天津大學 管理與經濟學部,天津 300072)

0 引 言

車輛的碳排放量與行駛速度密切相關,受早晚交通高峰的影響,城市內車輛平均行駛速度在不同時段變化明顯,形成時變網絡,從而影響車輛的碳排放量[1]。城市車輛配送路徑問題(vehicle routing problem,VRP)已經越來越多考慮了時變網絡影響[2,3]。但是還需考慮新到達訂單對調度方案造成的擾動[4]。目前已有學者對客戶訂單的調度周期等實際配送問題進行研究[5],并對隨機到達的訂單擾動進行重調度來獲取更優方案[6]。車輛選擇不同路徑,平均行駛速度及行駛距離會有差異,碳排放量也會受到顯著影響[7,8]。蟻群算法在路徑優化過程中存在效率較低、易早熟等缺點,可通過改進算法流程和改進信息素更新方式,來提高算法效率,避免陷入局部最優[9]。

本文將以碳排放量最小作為主要優化目標,根據隨機訂單到達的時間安排調度,研究時變網絡、柔性路徑下配送貨車動態重調度的問題。即RTVVRP-PF(rescheduling of time-varying vehicle routing problem with path flexible)。采用模擬退火與蟻群混合的算法求解模型,以避免陷入局部最優,提高算法收斂速度;引入碳排放因子,改進信息素濃度更新方式;使用自適應精英個體繁殖策略提高算法效率。

1 問題描述與假設條件

1.1 問題描述

本文以城市內貨物配送為研究對象,貨運車輛從一個配送中心出發,為多個不同時間窗口的客戶需求進行配送。車輛可以在多條柔性路徑中,選擇到達客戶節點的最佳路徑,并且在時變網絡下,車輛速度隨道路擁堵程度、交通高峰期動態變化,車輛速度的變化能夠直接影響車輛的碳排放量大小。

(1)動態重調度對解空間的影響

調度池中新訂單的增加,將減少優化方案中的變量約束,擴大全局優化的解空間,從而找到更優的可行解,如圖1所示。

圖1中,原調度方案,車輛配送訂單2時,部分行駛時間處于交通擁堵期,車輛的速度降到25 km/h,由于速度較慢,車輛在交通擁堵期將產生大量碳排放。

假如在原調度方案中的車輛還沒有出發時,加入新訂單,一起優化調度,新訂單路徑平均速度較快(35 km/h),部分服務時間在交通擁堵期內不產生碳排放量,并且訂單2由于重調度,避開了高峰,速度變為55 km/h,將減少大量碳排放。因此,以一定周期加入新訂單進行多次動態重調度,比一次性的靜態調度具有更大的解空間。

(2)時變網絡下柔性路徑對碳排放量的影響

在現實的城市配送道路網絡中,網絡節點之間往往有多條路徑連通,即存在柔性路徑。在早晚交通高峰期影響下,不同路徑(路段)上,車輛平均速度變化較大。車輛在時變網絡的不同時段,通過選擇碳排放量較小的路段組成的路徑,可以減少碳排放量。

RTVVRP-PF模型以動態重調度的方式為基礎,考慮了在車輛平均行駛速度受交通影響的時變網絡下,將一次性的靜態車輛調度分成多個階段,即隨著客戶訂單的不斷到達,在積累一定時間或一定數量的訂單后,及時安排重新調度,通過新舊訂單的混合來重構解空間,以使得車輛配送盡量避開交通擁堵期,提高車輛平均行駛速度來減少碳排放量。并且為了貼合實際工況,在模型中考慮柔性路徑、動態載重、時間窗口約束等因素,并在模型中這些因素整體優化,以達到在VRP中減少碳排放的目的。

1.2 符號說明

模型中的變量定義如下:

V:節點集合,且V=1表示配送原點,V=(2, 3, …,N) 表示客戶節點集合;

E:網絡圖中的弧集,i,j為任意兩個客戶節點,E={(i,j)|(i,j∈V且i≠j)};

r:路徑上的路段,r∈{1, 2, …,R};

k:重調度次數,k∈{1, 2, …,K};

m:配送車輛, ?m∈M;

g,G:車輛的實時載貨量和最大載重;

qj:節點j的貨物需求量;

tj:車輛到達客戶節點j的時刻;

tsj:節點j的服務時間;

tsm:車輛m的裝載時間;

1.3 條件假設

(1)配送原點擁有足夠多的載重量相同的車輛。

(2)每次調度中,每個客戶只接受一輛車輛配送服務。

(3)車輛到達客戶的時間具有時間窗口限制。

2 模型構建

2.1 碳排放量計算

車輛碳排放主要由車輛行駛過程中所消耗的燃油產生,車輛燃油消耗估算的方法分為實驗法(MEET)與模型法(CMEM)兩類。模型法多采用燃油消耗量為目標建立模型,來估算碳排放量。但是模型法受到燃油不能完全燃燒的影響,車輛在不同速度下,燃油消耗量與碳排放量之間關系并不是恒定的,會發生變化,因此會影響時變網絡下計算碳排放量的準確度[10]。MEET法應用范圍較廣[11],本文采用MEET方法,即通過對貨運車輛的尾氣排放污染進行大量實驗分析,直接給出不同類型空載的貨車在不同速度v(km/h)下的碳排放量函數θu(g/km)和系數a,并可根據配送車輛類型和性能,對系數進行調整[12],計算公式表示為

(1)

為適合在城市內范圍進行配送,本文選用載重量為4.5噸的貨車,根據文獻[12],a0,a1,…,a6的取值依次為110,0,0,0.000 375,8702,0,0。

貨運車輛的碳排放量與載貨量成正比關系,文獻[12]中不同類型車輛在不同速度下,車輛滿載的碳排放量θl(v)的計算公式為

(2)

根據貨運車輛的載重量,b0,b1,…,b5取值不同,裝載量4.5噸車輛的取值依次為1.27,-0.002 35,0,0,-1.33。裝載量4.5噸車輛滿載與空載相比較,不同速度下,碳排放量系數的變化如圖2所示。

可以看出,4.5噸貨車滿載下,速度v高于23 km/h后,速度越快,碳排放系數越低,單位碳排放量越少。

綜上,某貨運車輛m在速度v下,實時載貨量gm,其動態實時碳排放量θ可以表示為

(3)

2.2 服務時間

服務時間包含貨運車輛裝車、卸車和與客戶進行手續交接的時間,在服務時間內車輛不產生碳排放量。車輛裝載在配送原點完成裝載,因此其裝載時間tsm不計入車輛行駛調度模型。客戶j的服務時間tsj只需要考慮卸載時間tsj1和手續交接時間tsj2,卸載時間與客戶j的需求量qj和單位時間卸載量qt(噸/小時)相關,即:tsj1=qj/qt。 如果tsj1=0, 即沒有貨物卸載,則tsj=tsj2=0, 表示在沒有到達客戶節點時,服務時間為0。客戶j的服務時間可以表示為

(4)

2.3 模型建立

RTVVRP-PF重調度時變網絡下柔性路徑低碳城市配送模型中,選取3個優化目標:所有車輛一天配送過程產生的碳排放量總和C最小,車輛行駛的總里程D最小,車輛總行駛時間T最小。C為主要優化目標,車輛碳排放量和D密切相關,且D也為VRP主要優化目標,另外,由于存在配送時間窗口約束,且T與車輛司機、車輛日常管理成本密切相關,因此也非常重要。本文選取C、D、T作為模型的3個優化目標。

T、D、C優化目標函數分別為

(5)

(6)

(7)

目標函數約束條件S.T.

(8)

(9)

(10)

(11)

(12)

(13)

(14)

式(8)表示車輛最大負載限制;式(9)表示客戶時間窗口約束;式(10)表示車輛完成任一客戶e的配送任務后需要離開;式(11)表示車輛從配送原點出發一次,完成全部配送任務后需要返回原點;式(12)表示路徑路段的決策變量;式(13)表示選擇客戶節點的決策變量;式(14)表示重調度決策變量。

3 算法設計

VRP屬于旅行商問題,早被驗證為NP完全(non-deterministic polynomial complete)問題,難以使用數學解析的方法進行求解[13]。智能優化算法已經公認對此類NP完全問題有較好的優化效果。蟻群算法(ant colony optimization,ACO)和模擬退火算法(simulated annealing,SA)在VRP優化求解中效果較好[14,15]。ACO是一種群智能優化算法,具有記憶性,利用螞蟻的信息素濃度變化來選擇和構造可行解。ACO每次都需要對整個種群重新構造,因此效率較低,而且容易陷入停滯。SA以某種概率接受更新解,容易跳出局部最優和停滯,且控制簡單,編碼易實現,但是計算時間較長,尋優效果一般,優化效率不高。

本文將上述兩種算法的優點結合,提出混合蟻群算法SA-ACO,算法通過對SA中退火溫度的控制,實現快速收斂,并利用SA的Metropolis準則防止陷入局部最優和過早停滯現象。通過ACO中的多參數信息素濃度計算來更新種群,增加種群的記憶性,并設計了碳排放因子,以增強種群優化的方向性。SA-ACO的算法流程如圖3所示。

跳出內循環的判斷標準為達到設定的某一個循環次數的常數值Cr1,跳出外循環標準為設定的某一個循環次數常數值Cr2內最優解的適應度不再發生變化。

SA-ACO算法中按照Metropolis抽樣準則來接受個體si(t)是否更新為s′i(t), 避免早熟或停滯。接受概率函數Ω為

(15)

其中,Г為退火為溫度控制參數,其迭代計算見式(16),L為算法的迭代次數

(16)

3.1 算法編碼

SA-ACO算法采用實數方式編碼,每個可行解為一個兩行的矩陣,矩陣的第一行為配送客戶節點順序的信息,數字“1”表示配送原點,車輛從原點出發,完成配送任務后返回原點。可行解矩陣第二行為柔性路徑選擇信息,并且部分節點之間的路徑包含多個順序排列的路段,如圖4所示。

圖4第一行編碼中 [1,4,3,11,1] 表示一個配送車輛從原點“1”出發,順序到達客戶節點4、3、11進行配送服務,最后返回原點。第二行路徑選擇信息,即車輛在1、4節點之間選擇了第2條路徑,4、3節點之間選擇了第1條路徑,3、11之間選擇了第3條路徑,從客戶節點11返回原點選擇了第2條路徑。其它配送車輛行走信息依次類推,直到所有客戶配送服務均完成。

3.2 適應度計算

SA-ACO算法中,種群個體的適應度與優化目標相關,模型包含碳排放量C、總行駛路徑D和總行駛時間T這3個優化目標,存在Pareto解集。3個優化目標單位不同,且數量級相差較大,在個體適應度計算過程中需對優化目標的數量級進行調整,以便適合其權重系數。本文采用自適應的適應度函數f(s)來消除多個目標數量級的差別,f(s)表示為

(17)

其中,Cmax(S)、Dmax(S)、Tmax(S) 表示某一次迭代種群中,個體s的碳排放量、總行駛路程和總行駛時間的最大值,λ1、λ2和λ3分別表示三者權重系數,需要根據企業對三者的關注程度決定權重系數。

3.3 個體更新

本文提出自適應精英個體繁殖策略來更新個體,來保證能夠搜索到全部解空間,動態保留精英個體的優良基因,以提高算法進化效率。

SA-ACO的種群個體采用兩點局部更新的方式,如圖5所示。

在客戶配送節點信息編碼中,將每臺車輛出發的原點設為編碼1,找到所有的U個車輛配送原點“1”,隨機選取第u1和u2兩個配送原點,將兩個配送原點中間的編碼按照ACO中信息素濃度進行重新構造,兩個隨機點之外的其它部分的信息編碼直接復制,構成新的種群個體。其中這兩個點之間的長度u12的計算公式為

(18)

式(18)中int[ ]向表示下取整,fmax(S) 為這一代群體中的最大適應度值,并且如果u12>U-u1, 則u2為個體的最后一個配送原點,即u1之后的部分全部更新。

自適應精英個體繁殖策略,能夠快速搜索到全部解空間,按照信息素濃度構造個體局部,能夠使其更新具有方向性的同時,保留其它部分的優良信息。

為增加模型中的碳排放量優化目標效果,并兼顧總行駛里程和總行駛時間的優化目標,本文在構造ACO的信息素能見度和信息素增量更新時,設計了包含碳排放量因子的多因子計算方式。

(19)

其中,allowedj′表示還沒有配送的客戶節點集合。

(1)包含碳排放因子的多因子能見度更新

(20)

(2)包含碳排放的多因子信息素增量計算

(21)

其中,0<σ≤1,為信息素蒸發速率。

(22)

其中,ε1、ε2、ε3為權重系數。

4 仿真算例

4.1 算例數據

為驗證RTVVRP-PF模型與算法的有效性,盡量貼近現實城市配送狀況且不失一般性,本文依據我國城市的配送數據的統計分布,利用計算機生成符合統計分布的隨機算例數據。

城市內配送,客戶的需求量一般較小,許多城市限制配送貨車的載重量,一般不使用較大型貨車。由于5噸以下的異質貨車(車型不同、載重量不同的貨車),碳排放量差異不大。所以這里使用最大載重量G=4.5噸的同質貨車(車型一致的貨車)進行配送運輸的模擬計算。

表1 隨機節點的坐標、需求量(噸)和時間窗口

由于相同時刻,城市內各條道路的車流密度不同,車輛的平均行駛速度也不盡相同,有些城市的快速路和擁堵路段與其它道路的平均速度差異較大。隨機選取10%的路徑,按照較慢速度[20~40] km/h和較快速度[55~70] km/h,產生以5為步長的整數平均分布作為車輛行駛速度。兩節點間不同路段,速度差異較大的路徑,則按照均勻分布隨機產生速度和路段長度百分比,見表2。

表2 柔性路徑的速度/(km/h)和路段長度/%

利用交通出行網以及地圖軟件,統計整理一般城市各個時段道路的平均行駛速度,見表3。表2中的路徑為柔性路徑,所有柔性路徑的平均速度為60 km/h。按照表2中柔性路徑不同時段的平均速度與柔性路徑平均速度60 km/h的比值,來確定柔性路徑在不同時段的速度系數(見表3)。則柔性路徑在不同時段的速度為表2中的速度乘以不同時段下的柔性路徑速度系數。

柔性路徑下,兩節點之間多條路徑的長度不同,這里將第二、三條路徑按照第一條路徑的長度以分布函數按照一定比例進行放大,放大的比例服從U(1,1.3)平均分布,見表4。

表3 車輛不同時間段的平均行駛速度

表4 柔性路徑的第二、三條路徑長度/km

4.2 仿真結果

本文按照4.1節中的數據進行仿真運算。使用Matlab2019a版本進行建模優化,種群規模取500。參數設置直接影響算法求解的效率與效果,本文采用均勻設計試驗法來確定參數。算法參數中退火初始溫度Г0=1500,降溫系數0.99,內循環次數Cr1=1000,算法終止(跳出外循環)條件為Cr2=1000次迭代內,算法最優解適應度不再變化。算法的路徑選取概率函數中α和β的數值取值為1,信息素的蒸發率σ=0.58。優化目標權重的系數λ1、λ2和λ3取值均為1,即企業對3種優化目標的重視程度相等。

重調度過程中,需要統計訂單的到達時間,假設7點半之前,前一日未配送訂單和新增的訂單為客戶節點2~15的訂單,并且7點半之后,每1.5小時,新增5個客戶訂單,單位貨物的卸載時間qt=0.35噸/小時,每個客戶節點的手續交接時間tsj2=0.2小時,重調度后,揀貨和裝車的時間tsm=1,訂單到達時間見表5。

表5 訂單到達時間

分別采用動態重調度和靜態調度方式進行仿真對比運算。首先采用動態重調度方式,即按照表5的訂單截至時間進行4次調度,并且每次重調度,都將上一次沒有出發的貨車和訂單重新加入新的調度方案,調度完成后,經過1小時揀貨和裝車,貨車可出發。如果重調度時貨車已開始裝車,則該貨車和訂單不能加入新的調度方案。然后采用靜態調度方式,即只進行兩次調度,在7點半進行2~15客戶節點訂單調度,在12點進行16~30客戶節點訂單調度,并且第2次調度時,不加入第1次調度中沒有出發的貨車和訂單。分別采用兩種調度方式,各進行30次優化運算,取平均值,結果見表6。

表6 動態重調度與靜態調度結果對比

由表6可以看出,與靜態調度方式相比較,采用動態重調度方式,在碳排放、總行駛時間、總行駛里程和使用車輛數方面均有所降低,尤其碳排放量C(kg)降低了15.3%,這主要是因為重調度使得優化空間加大,更多的車輛可以避開交通擁堵期,從而使得碳排放量等指標能夠明顯減少。

以30次優化計算中的一個重調度方案為例,重調度優化方案見表7,優化后的碳排放與總行駛時間和里程見表8。

表7 重調度優化方案

由表7可以看出,在當日7點半時,對客戶節點第2~15的訂單進行調度,共安排4輛貨車,最早貨車8點半裝載貨物后出發。截至到當日上午9點時,接到客戶節點16~20的新訂單,這些訂單與第一次調度時尚未開始裝車的訂單(分別安排在9點半、10點半裝完車),即客戶節點2、4、5、6、9、12、15的訂單一起進行第2次重調度,然后在上午10點半與上午12點,以相同方式分別進行了第3次和第4次重調度,共使用貨車8輛。

圖6為車輛配送路徑,可以看出,部分路徑有交叉,主要原因是訂單配送是經過多次重調度產生的,以及優化目標中增加了減少碳排放目標。

多目標函數的Pareto解與各優化目標的權重系數相關,RTVVRP-PF的3個優化目標可以根據企業對優化目標的需求程度來確定。本文對3個優化目標給出參考范圍。創建Pareto非支配解集合,根據支配關系進行Pareto分級,將級別最高的解放入Pareto非支配集合中,若非支配解集規模超過設定值,則進行修剪。利用上述方法求解得到10個Pareto非支配解,見表8,并分別求出了3個優化目標的極小值分別為碳排放252.45 kg,總行駛時間33.28小時和總行駛里程631.36 km。

表9為這次調度的碳排放、總行駛時間和里程。并且,為驗證柔性路徑和動態負載對配送優化目標的影響,在這個調度方案中,分別統計了相同調度方案下,如果使用最短單一路徑和貨車空載、滿載的情況下,車輛的碳排放量、總行駛時間和里程變化。

由表9可以得出,如果不使用柔性路徑,而使用單一路徑,并且是最短路徑,優化后,雖然總行駛里程稍有降低,但是碳排放量增加了11.2%,總行駛時間增加了9.2%,因為在柔性路徑狀況下,車輛在道路擁堵時可以選擇速度更快、碳排放量更低的較長路徑,說明柔性路徑更加符合實際配送工況。由表9還可以得出,不考慮動態負載的影響,分別直接使用固定的貨車空載和滿載質量來計算碳排放量,與考慮動態負載相比,碳排放量分別相差7.5%和6.5%,表明動態負載對貨車碳排放量的影響是比較顯著的。

表8 Pareto非支配解集

表9 柔性路徑及動態負載對優化目標的影響

為驗證分路段對城市配送碳排放量的影響,計算表7優化方案中的3個包含分路段的路徑碳排放量。當分路段時,分別計算每個路段的碳排放量并求和,當不分路段時,按照整條路徑的平均速度計算碳排放量,計算結果見表10。可以看出,城市配送中,分路段對碳排放量影響較大,尤其是各路段速度差異較大時,而同一時間,城市道路擁堵程度往往差異較大,所以在城市配送調度優化過程中,不能忽略分路段對優化結果的影響。

表10 分路段對優化目標的影響

為驗證SA-ACO算法的效果,分別使用粒子群(particle swarm optimization,PSO)、SA和ACO算法,按照4.1節隨機算例數據,采用CPU3.3GHz,內存8 G,64位Windows10操作系統的計算機,SA中采用實數編碼,降溫系數取0.99,初始溫度1500,每個溫度迭代次數1000;PSO中慣性權重0.73,兩個加速因子均為1.21;ACO算法的參數取值與SA-ACO相同。各計算30次,結果取平均值。均采用均勻設計方法設置算法參數,結果見表11。

由表11可知,SA-ACO在碳排放量、總行駛時間、總行駛里程和使用運輸車輛數量方面實驗結果的數值均為最小。在計算耗時方面,由于SA-ACO混合蟻群算法在構造路徑方面計算量較大,相比較于PSO算法中尋找粒子方向的更新方式,SA-ACO計算耗時略高于PSO算法。

表11 算法效果比較

由圖7可知,與其它算法相比較,SA-ACO的收斂速度在100次迭代后,明顯優于其它算法。

為進一步分析SA-ACO算法效率,本文采用TSPLIB中的3個經典的算例進行優化分析,算法參數的取值與上述RTVVRP-PF模型相同,優化結果見表12。

表12 TSP模型算法效果比較

由表12可以看出,SA-ACO對TSP以及VRP問題的優化計算具有較好效果。

5 結束語

在城市配送RTVVRP-PF模型中,首先驗證了動態重調度和靜態調度相比較,對減少碳排放量的效果明顯。其次,提出SA-ACO算法,提高了算法效率和跳出局部最優能力;設計了自適應精英個體繁殖策略,提高了種群優良基因的個體數量;引入包含碳排放的多因子算子,增強信息素更新的方向性。最后,按照城市配送實際工況,隨機生成算例,通過對算例的運算結果進行分析,驗證了模型在動態重調度、柔性路徑、動態負載和分路段方面,對減少碳排放量效果顯著,同時驗證了SA-ACO算法的高效性。

猜你喜歡
優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
PEMFC流道的多目標優化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
圍繞“地、業、人”優化產業扶貧
今日農業(2020年16期)2020-12-14 15:04:59
事業單位中固定資產會計處理的優化
消費導刊(2018年8期)2018-05-25 13:20:08
4K HDR性能大幅度優化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 成人在线视频一区| 欧美中文字幕一区| 久久性妇女精品免费| 91色综合综合热五月激情| 精品一区二区三区自慰喷水| 五月婷婷精品| 欧美有码在线| 色噜噜中文网| 四虎永久在线精品国产免费| 日韩欧美网址| 超碰色了色| 国产超薄肉色丝袜网站| 国产打屁股免费区网站| 欧美一级在线看| 亚洲精品中文字幕午夜| 国产午夜福利在线小视频| av在线无码浏览| 色首页AV在线| 精品国产污污免费网站| 亚洲一区网站| 日韩午夜伦| 国产精品久久久久久久伊一| 国内精品91| 五月天在线网站| 亚洲综合亚洲国产尤物| 日韩精品一区二区深田咏美| 欧美亚洲一区二区三区导航 | 亚洲中文无码h在线观看| 97视频免费看| 国产精品福利在线观看无码卡| 色婷婷亚洲综合五月| 在线观看国产一区二区三区99| 97av视频在线观看| 丁香综合在线| 成人福利在线免费观看| 在线观看av永久| 精品日韩亚洲欧美高清a| 亚洲天堂网视频| 国产69囗曝护士吞精在线视频| 国产精品免费露脸视频| 欧美色视频在线| 国产丰满成熟女性性满足视频| 日韩精品高清自在线| 久草视频精品| 激情亚洲天堂| 免费视频在线2021入口| 色婷婷综合在线| 日韩欧美成人高清在线观看| 久久精品国产999大香线焦| 国产精品播放| 狠狠干综合| 欧美中文字幕无线码视频| 日本手机在线视频| 欧美一级夜夜爽www| 国产丝袜第一页| 性色一区| 九九久久精品免费观看| 99re66精品视频在线观看| 成人精品视频一区二区在线| 亚洲热线99精品视频| 国产成人综合久久| 欧美性色综合网| 国产真实乱了在线播放| 国产精选自拍| 成人福利在线看| 国产区精品高清在线观看| 韩国福利一区| 97se亚洲综合| 毛片在线播放网址| 暴力调教一区二区三区| 免费国产高清视频| 最新国产午夜精品视频成人| 欧美日韩北条麻妃一区二区| 巨熟乳波霸若妻中文观看免费| 国产原创第一页在线观看| 四虎影视库国产精品一区| 日韩经典精品无码一区二区| 久久99蜜桃精品久久久久小说| 少妇高潮惨叫久久久久久| 激情综合婷婷丁香五月尤物 | 国产在线日本| 91亚洲视频下载|