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

基于ADMM方法的新增列車條件下靈活的列車時刻表優化算法研究

2021-03-12 08:56:26高如虎牛惠民
鐵道學報 2021年2期
關鍵詞:優化方法

高如虎,牛惠民

(蘭州交通大學 交通運輸學院,甘肅 蘭州 730070)

在高速鐵路(以下簡稱高鐵)運輸組織中,當現有列車服務難以滿足突發增長的客流需求時,需要在原有運營列車時刻表的基礎上,規劃新增列車的時刻表。對于新增列車條件下的時刻表優化是“列車運行調整”問題,但實質也屬于列車時刻表優化范疇。列車時刻表問題(Train Timetabling Problem, TTP)是軌道交通領域中的經典問題。列車時刻表問題涉及數量眾多的車站和列車,模型中所含決策變量數目隨問題規模呈指數式增長,屬于典型的NP-hard問題,很難在多項式時間內找到一種有效算法對其精確求解。尤其對于不同速度等級列車共線運行的情況,還需要處理列車越行時間和位置選擇問題,使得問題更加復雜。如何設計可靠有效的求解方法,始終是列車時刻表問題的主要挑戰[1-2]。

構建整數規劃模型并利用商業優化軟件求解列車時刻表問題,是近年來興起的一種方法。文獻[3-5]構建列車時刻表的混合整數規劃模型,并利用優化軟件CPLEX求解。文獻[6]以總的乘客等待時間最少為優化目標,研究了時變客流需求條件下列車時刻表的優化問題,建立了帶有線性約束的非線性整數規劃模型并利用通用優化軟件GAMS求解模型。不同速度等級列車共線運行情況下,還需要處理列車的越行問題,該情形改變了列車的發車順序。因此,模型中除了列車到發時刻變量外,還需要表示列車發車順序的0-1變量,這使得決策變量和約束條件的數量急劇增大[7-8]。而利用商業優化軟件求解對數學模型要求苛刻,僅適合中小規模問題。

基于問題分解的求解方法,是解決列車時刻表問題的另一種有效思路,核心是將復雜的問題分解為容易求解的多個子問題,分別求解各子問題,然后通過集成和迭代,獲得問題的最優解。近些年,很多文獻利用拉格朗日松弛算法對復雜的列車時刻表問題進行分解并取得了顯著的求解效果。文獻[9]構建了基于時空圖的時刻表整數規劃模型,利用拉格朗日松弛算法將列車之間的耦合約束松弛,并分解為單列車在時空圖中的最短路徑子問題。文獻[10]在研究城市軌道交通的列車時刻表問題時,考慮了能力和資源約束,最后利用拉格朗日松弛算法求解。文獻[11]提出了基于累計變量的列車時刻表整數規劃模型并設計了拉格朗日松弛算法求解。文獻[12]研究了新增列車條件下的列車時刻表優化問題,并綜合考慮了車站進路,利用拉格朗日松弛算法分解原問題。ADMM在拉格朗日松弛方法的基礎上,通過引入增廣項可快速的獲得原問題的可行解。由于ADMM方法優越的分解性能,在大數據和機器學習領域[13]及分布式計算領域[14]被廣泛使用。最近,文獻[15]將ADMM方法應用到車輛路徑優化問題,并提出了ADMM在交通運輸組織優化領域的應用前景。文獻[16]構建了基于擴展時空網絡的周期性列車時刻表模型,并提出了基于ADMM方法的分解機制。

現有研究大多固定了列車的區間運行時間和列車停站方案或者為規劃列車指定相對較小的出發時間窗,又或者固定了列車的發車順序。本文以新增列車條件下的時刻表優化問題為切入點,創新點為:①為獲得更高質量和更符合實際的列車時刻表,提出了更加靈活的列車時刻表設置方法。具體地,列車可靈活的選擇區間運行時間、停站方案、停站時間、到發順序、越行時機,此外,為規劃列車賦予以小時為出發時段的更加靈活的出發時間范圍。這種處理方法更加靈活,也更符合實際。②這種靈活的列車時刻表優化架構需要設計一種更加有效和更加科學的方法來解決列車對時空資源占用選擇沖突問題。為此,本文提出了一種更加緊湊的基于時空弧段的不相容約束刻畫列車的耦合約束。③分別利用標準的拉格朗日松弛方法和ADMM方法解決新增列車條件下的靈活列車時刻表優化問題,并對兩種求解方法進行了比較,驗證了ADMM方法在求解該優化問題時的優越性。

1 問題分析及說明

1. 1 問題分析

本文研究雙線高速鐵路新增列車條件下的時刻表優化問題。按照運輸生產實際,鐵路2個方向列車相互獨立運行,選取其中一個方向進行研究。沿著該線路列車運行方向:U為車站集合;Ie為原有時刻表中列車集合;Ia為新增列車集合;[0,T]為全天運營時段,不失一般性,按照1 min為間隔對全天運營時段進行離散處理,以1 min的整數倍表示列車在車站的到發時刻和在區間的運行時間。

本文提出靈活的列車時刻表優化架構,主要有以下幾方面的考慮:首先,鐵路運營者往往按照小時客流需求推算運營列車數量,靈活架構下給定每個小時的列車數量與鐵路運營實際相契合。其次,列車在線路區間運行速度并不固定,這種區間運行時間的靈活設置與鐵路運營實際更加契合。然后,列車停站方案與列車時刻表關系密切,應該對二者進行綜合設計。最后,傳統的指定列車較小出發時刻范圍的方法,事實上限制了列車對出發時刻的選擇,并且相對固定了列車的出發順序,此外,列車的出發時刻范圍在實際中也很難確定。

事實上,在鐵路運營實踐中,這種靈活的列車時刻表優化架構相對于傳統的列車時刻表優化更加貼合實際。理論上,這種靈活的優化方法相當于傳統列車時刻表優化對列車出發時間窗、停站方案、停站時間以及列車發車順序的松弛處理,此種處理能夠獲得更高質量的列車時刻表,但也無疑增加了列車時刻表優化的復雜度,需要設計一種更加有效的方法對其求解。

1.2 符號說明

本文通過構建時空網絡描述列車時刻表問題,并據此建立基于時空弧段的0-1整數規劃模型。以下給出用于構建網絡和建立模型所需要的集合、索引、參數和變量。

索引:i、j為列車索引,i、j∈I;u、u′為車站索引,u、u′∈U;t、t′、τ、τ′為時刻索引,t、t′、τ、τ′∈[0,T];(u,t)、(u′,t′)為時空網絡節點索引,(u,t),(u′,t′)∈N;(u,t;u′,t′)為時空網絡中從節點(u,t)到(u′,t′)的弧段索引,(u,t;u′,t′)∈A。

決策變量

xi(u,t;u′,t′)=

2 時空網絡構建

由于列車時刻表固有的時空特性,構建時空網絡模型是解決列車時刻表問題最常用的建模方式。通過構建時空網絡,可將列車時刻表信息(如始發時刻信息、停站信息等)圖形化地顯示在網絡結構中。此外,這種表現形式還可以將時刻表問題中難以描述的約束條件,轉換為網絡中弧段或節點的制約關系(如可達性、不相容性、網絡流量守恒等),將優化目標轉換為列車占用網絡弧段的最小費用。G=(N,A)為時空網絡,對每一個列車,都對應一個子網絡。Gi=(Ni,Ai)為列車i對應的時空子網絡,其中,Ni為列車i可占用的時空節點集合,Ai為列車i可占用的時空弧段集合。下面以列車i為例說明時空網絡的構建過程,見圖1。

圖1 時空網絡圖

2.1 網絡節點

為保證時空網絡的完整性,對列車i分別引入虛擬起點σi和虛擬終點τi。對于每個列車i∈I=Ie∪Ia,頂點Ni集合包含虛擬起點σi、車站到達節點集合Ru、車站出發節點集合Wu和虛擬終點τi。值得說明的是,對于原有列車i∈Ie,經過各個車站對應的起始節點和到達節點集合可由該列車原始時刻表、始發時刻最大偏移量、最小最大停站時間、最小最大運行時間確定;對于新增列車i∈Ia,經過各個車站對應的起始節點和到達節點集合可由出發小時時段、最小最大停站時間、最小最大運行時間確定。列車i的時空節點集合Ni可表示為

Ni={σi,τi}∪{Woi,Woi+1,…,Wdi}∪{Roi,Roi+1,…,Rdi}

2.2 網路弧段

列車占用弧段的費用見表1,α,β分別為原有列車和新增列車的費用權重。原有列車和新增列車對應的優化目標不同,對于原有列車,期望新優化的列車時刻表與原有時刻表的偏離程度越小越好,包括在始發站發車時刻的偏離以及區間運行時間、車站停站時間的偏離。而對于新增列車,期望總的旅行時間最小,包括區間運行時間以及在車站的停站時間。具體地,對于原有列車i∈Ie,起始弧的費用表示選擇始發站的出發節點對應時刻與原列車時刻表始發時刻的偏差費用,運行弧和停站弧費用分別表征列車選擇該運行弧段和停站弧段對應的運行時間和停站時間與原有時間的偏差費用。而對于新增列車i∈Ia,可在出發小時時段任意時刻出發,起始弧的費用為0。對于列車占用運行弧和停站弧的費用分別表征列車在區間的運行時間和車站停站時間費用。列車占用終止弧的費用均設置為0。列車從虛擬起點到虛擬終點對時空網絡弧段的占用形成了一條時空路徑(如圖1紅色線條所示),該時空路徑反映了列車在各站的到發時刻信息。

表1 弧段費用

3 模型建立

3.1 優化目標

考慮到本文對列車區間運行時間、車站停站時間及新增列車出發時間窗靈活的設置,容許列車自由選擇停站、越行時機、區間運行時間。這些因素直接決定了列車的旅行時間,進而影響鐵路服務質量。因而,本文以新增列車的旅行時間最短為優化目標,同時我們希望原有列車的時刻表偏移量最小。通過時空網絡的構建,列車時刻表問題轉換為列車占用時空資源費用最小的網絡優化問題,目標函數為

(1)

3.2 約束條件

(1) 網絡流平衡約束

流平衡約束是網絡優化問題最基本的約束。通過時空網絡構建,在列車時刻表問題中,對于任意列車i,需要保證僅有一條時空路徑被列車選擇。流平衡約束為

(2)

(2) 列車停站約束

在原列車時刻表中,對于原有列車i∈Ie,若在u站停站,則在新優化列車時刻表中也必須在u站停站,該列車在其他車站可自由選擇停站。對于新增列車i∈Ia,根據途徑車站的類型,可分為兩類情況,其中一類為必停站,此類車站一般為始發終到站及具有大客流的車站;另一類為備選停站,列車經過此類車站時,可自由選擇是否停站及停站時間,此類車站一般為中間站及越行站。因此,列車停站約束可描述為

(3)

如果列車在某車站停站,則停站時間在最小和最大停站時間范圍內自由選擇。通過時空網絡的構建,停站時間約束已經在時空網絡停站弧中得到保證,這里再無需考慮。

(3) 安全間隔及越行約束

列車在時空網絡中除滿足上述的流平衡約束和列車停站約束外,還需要滿足列車安全間隔約束。此類約束是多列車之間的耦合關系,為了保證列車對于時空資源占用的排它性,本文引入不相容弧集合描述。

為保證列車運行安全,任意兩列車應滿足安全出發間隔和安全到達間隔約束。此外,不同速度等級列車需要滿足越行約束。事實上,列車安全到達間隔約束和列車越行約束均可轉換為列車在車站應滿足的出發間隔時間約束。為方便對列車安全間隔約束進行統一處理,本文引入不相容弧集合Φ(u,t;u′,t′),表示當某一列車占用時空網絡運行弧段(u,t;u′,t′)∈Atravel時,其他列車不能同時占用時空網絡中弧段的集合,即違反安全間隔的弧段。

給定某一時空弧段,則該弧段對應的不相容弧段集合也隨之確定。具體地,對于時空運行弧段a=(u,t;u+1,t′)∈Atravel和b=(u,τ;u+1,τ′)∈Atravel定義:

ha,b=max{hdep,(t′-t)-(τ′-τ)+harr}

(4)

表示時空弧段a和b若同時被列車占用,應滿足的安全間隔。

對于時空弧段a,如圖2(a)所示,若t≤τ

對于時空弧段a,如圖2(b)所示,若t-hb,a<τ

圖2 不相容弧示意

綜上所示,弧段a的不相容弧集合Φ(a)為

Φ(a)={b=(u,τ;u′,τ′)∈Atravel|t-hb,a<

τ

(5)

對于任意的時空運行弧段a,在不相容弧集合Φ(a)內,當有列車占用時空網絡弧段a時,其他任意列車不能占用不相容弧集合Φ(a)內的其他弧段。不相容約束可描述為

?i∈Ia∈Atravel

(6)

下面分兩種情況討論不相容約束:

(4) 變量域約束

xi(u,t;u′,t′)∈{0,1} ?i∈I(u,t;u′,t′)∈Ai

(7)

4 基于拉格朗日松弛的分解

本文所建立模型中的約束條件式(6)為多列車耦合約束,下面利用拉格朗日松弛算法對列車時刻表模型進行松弛和分解。具體地,引入非負的拉格朗日乘子λi(a)≥0,?i∈I,a∈Atravel,將模型中的耦合約束條件松弛,形成拉格朗日對偶問題。拉格朗日函數為

(8)

進一步地,拉格朗日函數可被整理為

(9)

其中

(10)

移除式(9)中的常數項,原問題被分解為如下的單列車子問題

(11)

s.t. 式(2)、式(3)、式(7)。

子問題為單列車在時空網絡中的最短路徑問題,可利用最短路徑算法快速求解。在拉格朗日松弛方法中,通過迭代求解子問題,并更新拉格朗日乘子可獲得原問題的下界。

在迭代求解時,使用次梯度法更新拉格朗日乘子。

(12)

ηq=1/(q+1)

(13)

拉格朗日松弛算法如下所示:

Step6終止條件。如果ggap≤ε(最小gap值) 或者q>Q(最大迭代次數),則算法終止;否則,令q=q+1,返回Step2。

在求解子問題時,如式(10)所示,列車的最短路徑費用取決于時空弧段的旅行時間費用ci(a)和拉格朗日乘子λi(a)。在每次迭代時,僅通過列車占用弧段的拉格朗日乘子值調節列車的路徑選擇。事實上,在拉格朗日松弛方法迭代過程中,當所有列車時空路徑選擇完成后,才統一更新各弧段的拉格朗日乘子值,相當于在一次迭代時同一小時內的同等級列車占用弧段的拉格朗日乘子也相等。這就導致了在每次迭代時,同一小時內的同等級列車會選擇同樣的時空最短路徑。這種解的對稱性問題嚴重阻礙算法的效率。

5 基于ADMM的求解方法

ADMM方法適用于可分結構的凸優化問題,可將原問題等價的分解成易求解的子問題,交替求解子問題,以獲得原問題的最優解或者可行解。ADMM方法是在拉格朗日松弛方法的基礎上引入二次懲罰項,構造增廣拉格朗日函數,將原問題分解成若干個容易求解的子問題;在次梯度迭代的過程中,按照一定次序對不同的變量進行交替迭代求解,最終得到最優或可行解。

5.1 基于ADMM的松弛、增廣和分解

根據ADMM原理,首先引入松弛實數變量di(a)∈[0,1],將耦合約束不等式(6)轉化為等式

(14)

對于本文所建立數學模型,在拉格朗日松弛函數的基礎上,引入二次懲罰項ρ,構造增廣拉格朗日函數

(15)

增廣拉格朗日函數中含有二次項,使得問題失去了可分解性。下面對增廣拉格朗日函數進行線性化處理,令

(16)

[di(a)+pi(a)-1]2}

(17)

顯然,由于存在松弛變量di(a),上式中還存在二次項,下面確定松弛變量di(a)的取值。將上式中除松弛變量di(a)以外的其他變量看成為參數,并將與松弛變量di(a)相關項組合為如下的二次函數形式

Ldi(a)=di(a)2+2×di(a)×[pi(a)-1+xi(a)+

(18)

式中:Q為與松弛變量di(a)的無關項。由于松弛變量di(a)為連續變量,因此,函數Ldi(a)取最小值的點為

(19)

下面分兩種情況進行討論松弛變量di(a)的取值:

(1) 當pi(a)≥1時,則di*(a)≤0,并且在[0,1]范圍內,函數Ldi(a)為遞增函數。故最優的d*(a)=0,帶入公式(18),二次項可轉化為

xi(a)+2×[pi(a)-1]2

(20)

綜合兩種情況,增廣拉格朗日問題可進一步分解為如下的單列車的子問題

(21)

s.t. 約束式(2)、式(3)、式(7)。

其中,

(22)

5.2 基于ADMM方法的求解過程

本文研究的列車時刻表問題含有關于多個列車的多個決策變量,將ADMM方法擴展為多塊(multi-block)的情況。值得說明的是,標準的two-block ADMM 方法的收斂性已經被證明[17]。而擴展的multi-block ADMM 方法被通過反例證明并不必然收斂[18]。盡管沒有收斂性的保證,但由于multi-block ADMM 方法的可操作性及高效性,在一些實際的優化問題中卻被廣泛應用[19]。

通過對原問題進行松弛、增廣、線性化和分解求解方法的設計,利用ADMM方法求解列車時刻表的具體流程如下所示。

Step2基于優先權迭代順序產生ADMM解。

For 每個列車i∈I

最短路徑算法求解i列車對應的ADMM子問題Lρ,i;

固定列車i對應子問題的解;

Endfor

Step6終止條件。如果ggap≤ε(最小gap值) 或這q>Q(最大的迭代次數),則算法終止并令上界解作為最后輸出結果;否則,令q=q+1,返回Step2。

值得說明的是,算法中Step2拉格朗日乘子和二次懲罰項的更新參考了文獻[15]中的方法。算法中Step3的啟發式可行算法與拉格朗日松弛算法中的啟發式可行算法相同,可利用文獻[11]中提出的優先權算法求解。

6 算例

以武廣高鐵線路為例,此線路包含18個車站,為方便描述,車站依次編號為1~18。在此線路上運行兩種速度等級的列車,不同速度等級列車在區間的運行時間見表2。高等級列車最小/最大停站時間分別為2、5 min,低等級列車最小/最大停站時間分別為2、10 min。線路上在全天規劃時段06:00—24:00內,原有列車135列,計劃新增列車20列。為便于處理,假定原有列車均為高等級列車,新增列車均為低等級列車,假定列車在區間的運行時間為常數。原有列車的始發站出發時間最大偏移量為10 min,關于新增列車的相關運營參數見表3。本文所提出的算法在Visual Studio 2013平臺上使用C++語言實現,算法運行環境為1臺CPU Intel Core i5-4 590 s 3.00 GHz, 4 GB內存的個人計算機。

表2 列車在線路區間的運行時間

表3 新增列車相關運營參數

分別利用標準的LR和ADMM算法求解,共迭代100次,兩種算法對應的上下界迭代過程見圖3,對應的優化ggap值變化情況見圖4。對于ADMM方法,最優的上界為3 774,最優下界為3 254,最優ggap值為13.78%,計算時間為302.74 s;對于LR算法,最優的上界為4 097,最優下界為3 287.33,最優ggap值為19.76%,計算時間為195.83 s。從兩種算法的迭代過程及計算結果中不難看出,盡管ADMM方法的計算時間略大于LR方法,但求得的最優上界和ggap值更優。此外,ADMM方法收斂的速度比LR算法更快。因此,對于本文優化的新增列車條件下的靈活列車時刻表優化問題,ADMM相比于LR具有更好的計算性能。

圖3 LR和ADMM算法的上下界迭代過程

圖4 LR和ADMM算法的ggap迭代過程

利用ADMM方法優化的列車運行圖見圖5。限于篇幅,對應的列車時刻表(整數形式,如0對應06:00,依次類推)上傳至https:∥github.com/ruhugao/Ahu。從運行圖中可以看出,部分新增的低等級列車可被原有的高等級列車越行,如新增的從武漢站開往廣州南站的列車1分別在郴州西和韶關站被原有高等級列車越行;部分原有列車改變了原有的列車時刻表,如原有列車4在其始發站武漢站的出發時刻為07:33,優化后的出發時刻調整為07:37。部分新增列車除了在計劃的列車車站停站外,在其他車站也可停站,可以在最小和最大停站時間之間靈活的選擇停站時間,如新增列車1在計劃之外的郴州站停站4 min。本文提出的這種靈活架構下,利用兩種分解算法求解時,可以迅速的得到可行解。

圖5 優化的列車運行圖

7 結論

本文針對新增列車條件下的高速鐵路列車時刻表進行優化。為了獲得更高質量的列車時刻表,提出了一種靈活的列車時刻表優化架構。通過構建時空網絡的手段,建立了基于時空弧段的0-1整數規劃模型,提出了一種更加緊湊的不相容約束用以刻畫列車時刻表優化問題的列車安全間隔約束。分別利用標準的拉格朗日松弛方法和交替方向乘子法對原問題進行分解和求解,并以武廣高速鐵路線路為例進行了實例分析。研究結果表明:

(1)相對于傳統的列車時刻表優化問題,在靈活的列車時刻表優化架構下,更能夠快速的獲得可行的列車時刻表。

(2)本文提出的更加緊湊的基于時空弧段的不相容約束能夠更加清晰的刻畫列車之間的耦合關系,同時這種不相容約束在ADMM方法對原問題分解過程中也更加方便。

(3)通過實例驗證,利用標準的拉格朗日松弛方法求解的最優ggap值為19.76%,而利用ADMM方法求得的最優ggap值為13.78%。兩種方法相比,盡管ADMM方法的計算時間略大于LR方法,但求得的最優上界和最優ggap值更優。此外,ADMM方法收斂的速度比LR算法更快。

猜你喜歡
優化方法
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 亚洲一区二区精品无码久久久| 一级毛片免费播放视频| 久久免费成人| 天天色天天综合网| 欧美日韩精品在线播放| 天天干天天色综合网| 日韩精品专区免费无码aⅴ | 97亚洲色综久久精品| 亚洲精品福利视频| 欧美亚洲国产精品久久蜜芽| 久久精品最新免费国产成人| 色吊丝av中文字幕| 国内精品视频在线| 国产女人爽到高潮的免费视频| 国产丝袜丝视频在线观看| 欧美中出一区二区| 久久国产av麻豆| 国产国语一级毛片在线视频| 日韩在线视频网| 精品视频第一页| 国产成人夜色91| 国产精品成| 国产视频自拍一区| 婷婷五月在线| 天堂在线www网亚洲| 国产a v无码专区亚洲av| 手机永久AV在线播放| 亚洲无码视频一区二区三区 | 亚洲天堂区| 伦伦影院精品一区| 无码网站免费观看| 精品国产免费人成在线观看| 青草精品视频| 国产成人精品视频一区二区电影| 国内老司机精品视频在线播出| av在线手机播放| 97se亚洲综合在线天天| 99精品视频九九精品| 亚洲日韩第九十九页| 伊人久久久久久久久久| 欧美一级高清片久久99| 亚洲精品国产成人7777| 91精品啪在线观看国产| 爱做久久久久久| 国产网站免费看| 在线亚洲精品自拍| 天天摸天天操免费播放小视频| 欧美日韩高清| 国产99在线观看| 2021最新国产精品网站| 午夜国产大片免费观看| 欧洲熟妇精品视频| 亚洲香蕉久久| 一本一道波多野结衣一区二区 | 欧洲日本亚洲中文字幕| 色网站在线视频| 国产精品中文免费福利| 69精品在线观看| 亚洲成A人V欧美综合| 精品精品国产高清A毛片| 国产美女主播一级成人毛片| 精品视频一区二区三区在线播 | 五月天久久婷婷| yy6080理论大片一级久久| 色老二精品视频在线观看| 亚洲免费毛片| 亚洲无码免费黄色网址| 五月婷婷综合网| 在线观看欧美国产| 伊人久久青草青青综合| 91色在线观看| 亚洲精品欧美日韩在线| 99在线视频免费| 国产va欧美va在线观看| 国产成人精品午夜视频'| 国产精品污视频| 亚洲福利片无码最新在线播放| 91偷拍一区| 中国国语毛片免费观看视频| 一级毛片在线播放免费观看| 欧美翘臀一区二区三区| 国产精品冒白浆免费视频|