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

考慮梯級水電耦合的水火電系統機組檢修計劃優化分解方法

2021-12-06 12:04:58代江田年杰姜有泉鄭志佳劉明波謝敏
電力大數據 2021年8期

代江,田年杰,姜有泉,鄭志佳,劉明波,謝敏

(1.貴州電網有限責任公司電力調度控制中心,貴州 貴陽 550000;2.華南理工大學電力學院,廣東 廣州 510640)

機組檢修計劃(Generator Maintenance Schedule, GMS)優化是指在滿足系統運行和機組特性等約束條件的前提下,確定一定調度周期內機組檢修任務開始時間的優化過程[1-2]。根據不同的優化目標,制定得出的檢修計劃有不同的效果,如提高系統可靠性[3-6]、保證經濟效益[7-9]等。而在清潔能源快速發展的背景下,考慮清潔能源的消納亦成為一種趨勢[10],如文獻[11]以最小化棄風量為目標,文獻[12]以棄水量最小為目標優化得出檢修計劃,文獻[13-14]綜合考慮了新能源的棄置量等。

我國西南地區水資源豐富,含有眾多的梯級水電站,梯級水電站之間存在復雜的耦合關系。因此,在建立水火電機組檢修計劃優化模型時,需要對梯級水電的耦合關系進行合理的描述。文獻[15]考慮梯級水電站上下游耦合特性對檢修計劃的影響,建立了梯級水電站的中長期調度和檢修計劃雙層優化模型。文獻[16]以發電收益最大化為目標,建立了梯級水電系統機組檢修計劃優化的混合整數線性規劃(Mixed Integer Linear Programming, MILP)模型,文中結合動態規劃、凸包逼近與窮舉法對水電站的出力進行了線性化。文獻[17]假設水電站出力是發電流量的一次函數,建立了考慮梯級耦合的水火電調度與機組檢修計劃協調優化MILP模型。

綜上,考慮梯級水電耦合的水火電機組檢修問題在數學上表現為一個含復雜約束的MILP模型,求解速度有待提高。采用分解法降低問題規模以簡化求解是一條有效途徑。分解法的基本思想是將原問題分為多個子問題,引入拉格朗日乘子將子問題之間的耦合約束增廣到目標函數中,從而將原MILP問題的求解轉化為各個子問題的求解。基于模型分解的拉格朗日乘子法的經典算法是拉格朗日松弛(Lagrange Relaxation, LR)法[18-19]。LR主要的缺點是:1)每一次迭代都需要求解所有的子問題以更新乘子,造成拉格朗日乘子出現振蕩現象;2)收斂速度較慢,在耦合約束多的情況下甚至難以收斂。增廣拉格朗日松弛法(Augmented Language Relaxation, ALR)[20]在目標函數中添加了二次懲罰項以加速收斂,但添加的二次項使得問題不可分。交替方向乘子法(Alternative Direction Multipliers Method, ADMM)[21-22]通過在求解時固定某一子問題的變量解決了這一困難,但無論是ALR還是ADMM,問題的形式都表現為混合整數二次規劃問題,其求解效率遠低于求解MILP的效率。

近年來,通過代理次梯度法更新拉格朗日乘子的研究日趨成熟[23]。該方法的優勢是在每一次迭代中只需求解部分子問題以獲得代理次梯度,而不需要求解所有子問題,從而節省計算時間。使用該方法更新拉格朗日乘子,能夠保證乘子的變化方向與下降方向構成銳角,即向收斂方向逼近。代理絕對值拉格朗日松弛法(Surrogate Absolute Value Lagrange Relaxation, SAVLR)[24]便是基于以上思想提出,該方法在代理次梯度法的基礎上,引入等價于二次懲罰項的線性懲罰項,改善收斂性能,線性懲罰項使得問題仍然保持線性性質,求解速度更快。

本文以最小化火電機組的運行成本、水電機組棄水成本和發電機組檢修意愿被調整的成本之和為目標,在滿足整個系統、水火電機組、水庫及輸電線路等的特性和約束條件下,建立考慮梯級水電耦合的水火電機組檢修優化MILP模型。再將原問題分解為火電和水電兩個MILP子問題,采用SAVLR法求解。并通過兩個水火電系統算例驗證了算法的有效性。

1 水火電機組檢修計劃優化模型

水火電系統的機組檢修計劃優化流程如圖1所示。由發電廠申報信息,電力調度中心根據負荷和天然來水的預測,按設定目標進行優化,制定合理的檢修開始時間,但不修改檢修時長,最后發電廠按電力調度中心下發的檢修開始時間安排檢修。

圖1 水火電系統發電機組檢修計劃優化的流程圖Fig.1 Flow chartof optimizing generator maintenance schedule in a hydrothermal power system

1.1 目標函數

檢修計劃優化的目標是最小化火電機組的運行成本、水電機組棄水成本和發電機組檢修計劃調整成本之和為目標,即:

(1)

目標函數包含兩項,第一項為火電機組運行成本與棄水成本之和。水電機組在運行時不會消耗燃料,發電成本為零,此處僅考慮火電機組的運行成本,其中Pi,t為第i個火電機組在t時刻的出力,ai、bi和ci分別是對應第i個火電機組成本二次函數的二次項、一次項與常系數;棄水成本即水電站計劃外排水的成本,其中ρ為電價,βh為流量—功率轉換系數,表示水電站h將流量轉化為出力的能力,qh,t為棄水量,表示水電站h當其庫容越限時需要排放的水量,該水量不進入水電機組發電,因此會產生相應的損失。

第二項為發電機組檢修計劃的調整成本。y為檢修計劃被調整的標志變量,當火電/水電機組i/j申報的檢修時間被調整時,yi/yj為1,否則為0;w為懲罰系數,取決于電力調度中心對發電廠滿意度的傾向程度。

模型中,T表示檢修時段集合,對于年度機組檢修計劃,本文以天為單位對機組檢修計劃進行安排,那么T=[1,2……365];I為火電機組集合;J為水電機組集合;H為水電站集合。

1.2 約束條件

約束條件包括功率平衡約束、火電機組與水電機組出力上下限約束、線路傳輸容量約束、旋轉備用約束、水電站的運行約束、梯級水電站上下游間的水量耦合約束、機組檢修的時長約束和連續性約束等。可以概括為三類:系統約束、水文約束及檢修約束,具體描述如下。

1.2.1 系統約束

1)功率平衡約束

整個系統的功率平衡約束可描述為:

(2)

其中,Pj,t表示在t時刻水電機組j的出力,J為水電機組數;Dt為系統在t時刻的負荷需求。

2)線路傳輸容量約束

輸電線路傳輸的功率都不應超過其上限,即:

(3a)

(3b)

其中,αi,l與αj,l為火電機組i與水電機組j對線路l的功率轉移因子;αn,l為負荷節點n對線路l的轉功率移因子;Dn,t為節點n在t時刻的負荷需求;N為負荷節點集合;Pmaxl,t為線路l的傳輸容量上限。

3)機組出力約束

(1-xi,t)Pi,min≤Pi,t≤(1-xi,t)Pi,max

(4a)

(1-xj,t)Pj,min≤Pj,t≤(1-xj,t)Pj,max

(4b)

其中,xi,t和xj,t為0-1變量,分別表示火電機組和水電機組的檢修狀態,為1時機組檢修,為0時機組不檢修,本文假設機組不檢修時處于啟動狀態;Pi,min與Pi,max表示火電機組的最小和最大出力;Pj,min與Pj,max表示水電機組的最小和最大出力。

4)旋轉備用約束

通過預留足夠的系統旋轉備用容量可以應對負荷預測誤差及事故帶來的影響,因此,添加旋轉備用約束可以保證系統運行的可靠性,即:

(5a)

(5b)

ri+rj=1+r

(5c)

其中,r為備用率。本文中r取為0.05。

1.2.2 水文約束

水文約束是針對水電站的約束,以下約束中,式(6)-(10)是對單個水電站內部變量的約束,式(11)為對上下游梯級水電站變量之間的耦合約束。

1)發電流量約束

0≤uh,t≤uh,max

(6)

其中,uh,t為第h個水電站在t時刻用于發電的水量,稱為發電流量;uh,max為第h個水庫的最大發電流量。

2)棄水量約束

0≤qh,t≤qh,max

(7)

其中,qh,t為水電站h在t時刻不通過水電機組排放的水量,該水量無法用于發電,稱為棄水量;qh,max為其上限。

3)水庫容積約束

vh,min≤vh,t≤vh,max

(8)

其中,vh,t為水電站h在t時刻水庫的庫容;vh,max和vh,min為水庫庫容的上下限。

4)水電站出力表達式

Ph,t=βhuh,t

(9)

其中,βh為水電站h的流量—功率轉換系數,表示將流量轉換為出力的能力,該值可利用水電站的歷史出力數據和歷史發電流量數據擬合得出。

5)水電站出力與機組出力之間的關系

水電站h的出力,等于該水電站內所有水電機組出力之和,如式(10)。

(10)

6)梯級水電站上下游水量耦合約束

上下游水電站之間的水量存在相互影響是梯級水電站一個重要的特點,在建模時需要進行適當考慮。通常在考慮上下游水量時,需要考慮水流滯時的影響,從而使得該式成為復雜的非線性關系式,但在本文中,由于以天為單位,故可以簡化對水流滯時的考慮,從而有式(11)。

vh,t-vh,t-1=s(fh,t-uh,t-qh,t+uh-1,t+qh-1,t

(11)

其中,s為將流量轉化為體積的系數;fh,t為水電站h處的天然來水,包含降雨和融雪等。

1.2.3 檢修約束

1)檢修時間約束

(12a)

(12b)

其中,di/dj表示火電/水電機組i/j檢修任務的持續時間,所有機組都必須完成規定時長的檢修任務。

2)檢修連續性約束

對于機組來說,檢修任務一旦開始就不能中斷,必須確保檢修任務的連續性,即:

xi,t-xi,t-1≤xi,t+di-1

(13a)

xj,t-xj,t-1≤xj,t+dj-1

(13b)

3)檢修調整變量約束

xi,ei-xi,ei-1≥1-2yi

(14a)

xj,ej-xj,ej-1≥1-2yj

(14b)

其中,ei/ej表示火電/水電機組i/j申報的檢修開始時間;yi/yj為0-1變量,表示火電/水電機組的檢修調整變量。

以火電機組為例,xi,t與yi,t之間的關系可用表1的真值表表示,若優化得出的檢修開始時間不同于申報的時間,則取為1,否則為0。

表1 x與y關系真值表Tab.1 Truth table of x and y

1.3 目標函數的線性化

對于目標函數中的二次項,本小節采用文獻[25]中的方法進行線性化,從而使原始模型轉化為混合整數線性規劃模型。線性化思路如圖2所示。

如圖2所示,已知火電機組i的最大出力Pi,max和最小出力Pi,min及其發電成本的系數ai、bi和ci,將其發電成本曲線分為兩段,其中間值為Pi,mid,分別計算出Pi,min、Pi,mid和Pi,max對應的發電成本Ci,min、Ci,mid和Ci,max后,將Ci,mid與Ci,min和Ci,mid分別連接,得到兩條直線(圖中紅色和藍色點虛線),設它們的斜率分別為Ki,1和Ki,2,截距分別為Mi,1和Mi,2,由此可得到發電成本Ci,t線性化后的表達式為:

Ci,t≥Mi,1xi,t+Ki,1Pi,t

(15a)

Ci,t≥Mi,2xi,t+Ki,2Pi,t

(15b)

圖2 發電成本曲線分段線性化Fig.2 Piecewise linearization of cost curve

通過添加如式(15)的約束,發電成本Ci,t一直在兩條直線上方,通過最小化發電成本,可以使發電成本與某條直線上的點重疊,達到將其線性化的目的。若要得到更加精確的逼近效果,可以通過將曲線分為更多段來完成。

(16)

線性化后,使用Ci,t替換(1)中的二次項表達式,可得如式(16)所示的目標函數。至此本文建立了一個水火電機組年度檢修計劃的MILP模型。

2 水火電機組檢修優化分解算法

2.1 模型分解

在求解大規模MILP問題時,可以采用分解的思想,將原問題分解為多個子問題求解,縮小問題規模,從而降低求解難度。對于本文建立的水火電機組檢修優化MILP模型,可以十分自然地將其分解為水電和火電兩個子問題,即:

火電子問題:

(17a)

s.t. (4a)(5a)(12a)(13a)(14a)(15)

(17b)

水電子問題:

(18a)

s.t. (4b)(5b)(6)-(11)(12b)(13b)(14b)

(18b)

火電子問題與水電子問題之間存在耦合約束如下:

(2)(3a)(5c)

(19)

為簡化表達,將火電子問題中的決策變量設為Xa,將水電子問題中的決策變量設為Xb。從而可以將火電子問題采用緊湊形式表示為:

minF(Xa)

(20a)

s.t.Xa∈Ωa

(20b)

類似地,水電子問題可以采用緊湊形式表示為:

minF(Xb)

(21a)

s.t.Xb∈Ωb

(21b)

火電子問題與水電子問題之間的耦合約束采用緊湊形式表示為:

g(Xa,Xb)=0

(22)

其中,式(20a)和式(21a)分別表示火電和水電子問題的目標函數;式(20b)與式(21b)為兩個子問題各自內部的約束,僅對屬于該子問題的決策變量起約束作用;式(22)為兩個子問題之間的耦合約束。

對于水火電機組檢修優化問題,耦合約束共有三組,分別是系統功率平衡約束(2)、線路傳輸容量約束(3a)和備用約束(5c)。對于不等式約束(3a),需要引入松弛變量將其改寫為等式約束:

(23a)

(23b)

其中,zl,t和ul,t為線路傳輸容量約束的松弛變量,恒大于零。因此,耦合約束(2)、(5c)和(23)就可以寫成式(22)表示的等式約束。

2.2 代理絕對值拉格朗日松弛法的應用

通過分解形成的火電子問題和水電子問題,無論在規模上還是在復雜程度上都遠小于原問題,易于求解。本節采用代理絕對值拉格朗日松弛法求解。

對于火電子問題(20),將耦合約束(22)松弛到(20a)中的同時,引入絕對值懲罰項(等價于二次懲罰項)改善收斂性能,則原優化問題(20)轉化為:

(24a)

s.t.Xa∈Ωa

(24b)

其中,λ為拉格朗日乘子;γ為恒正罰因子。

在目標函數中引入絕對值懲罰項可以加速收斂,但同樣也會使目標函數呈現非線性特征。為解決這一問題,對每一個絕對值引入一個連續變量及兩個約束將其線性化,具體形式如下:

(25a)

-om≤gm(Xa,Xb)≤om,m=1……M

(25b)

s.t.Xa∈Ωa

(25c)

其中,om為連續變量;M為耦合約束的數量。

以式(25)方式形成的子問題,在決策變量違反耦合約束時能夠給出足夠的懲罰,保證了收斂性能;同時也保持線性的性質,從而保證了求解的效率。需要注意的是,在求解火電子問題(25)時,應該固定水電子問題變量Xb;同理,在求解水電子問題時,也應該固定火電子問題變量Xa;而松弛變量在兩個子問題的求解中均作為變量進行求解。

在傳統的拉格朗日乘子法中,每一次迭代都需要求解所有子問題才能更新拉格朗日乘子。與之不同的是,本文使用的方法使用了代理次梯度的方式,在滿足代理最優條件的情況下只需要部分子問題的解就可以更新拉格朗日乘子,并保證更新方向與梯度下降方向成銳角,拉格朗日乘子的更新流程如表2所示。

表2 基于代理次梯度法的拉格朗日乘子更新流程Tab.2 Flow chart of updating Language multipliers based on surrogate sub-gradient method

對于火電子問題,其代理最優條件如下:

(26)

(27)

若求解某子問題后滿足代理最優條件,則使用該次迭代中子問題的值更新拉格朗日乘子,拉格朗日乘子更新公式如下:

(28a)

(28b)

(28c)

(28d)

(29)

在本文中,收斂判據ε取10-6。

2.3 算法流程

在求解考慮梯級水電耦合的水火電系統機組檢修計劃模型時,首先給變量及參數賦初值,并將問題分解為水火電兩個子問題,再使用SAVLR法迭代求解。SAVLR法的基本流程與拉格朗日松弛法的類似,均在每次迭代后,判斷耦合約束是否滿足收斂條件。不同的是在每一次迭代中,SALVR法依次求解2個子問題,若滿足代理最優條件(26)則直接利用式(28a)更新拉格朗日乘子;若在該次迭代中完成所有子問題的求解后仍不滿足(26),則需要減小罰系數重新開始該次迭代。分解算法流程如圖3所示。

圖3 分解算法流程Fig.3 Flow chart of decomposition algorithm

3 算例分析

本節采用3節點6機系統和某省級實際1151節點86機水火電系統為算例驗證基于SAVLR法的水火電機組檢修優化分解算法的有效性。

3.1 3節點6機水火電系統

在該系統中,節點2接入一個兩級的梯級水電站,節點1和3各自接入兩臺火電機組,該系統發電機特性與水庫特性取自某省級實際系統,火電運行成本分兩段進行線性化。

算例分別采用SAVLR和ADMM兩種方法求解。通過求可行解的方式賦初值,并將分解算法的結果與使用CPLEX求解器直接求解的結果進行對比,如表3所示。

表3 3節點6機系統中采用SAVLR和ADMM獲得的結果對比Tab.3 Comparison of results obtained through SAVLR and ADMM in 3-bus and 6-unit system

由表3可以看到,使用SAVLR法,約束被違反時會受到相應的懲罰,且懲罰項為線性形式,使得求解子問題時需要的計算時間較短,最終經過28.659秒,迭代7次后得到的目標函數值為986230.43,與CPLEX直接求解的結果相同。相比之下,ADMM盡管同樣對違反約束做出懲罰,但由于問題為混合整數二次規劃模型,同時因為每一次迭代中都需要求解所有的子問題,因此需要較長的計算時間,在244秒后經過13次迭代得到的目標函數值為989850.478。

3.2 某省級實際1151節點86機水火電系統

該系統包含621條支路、1151個節點(其中包含202個接負荷的節點)及41臺火電機組與45臺水電機組。其中擬安排檢修任務的有33臺火電機組和38臺水電機組。該系統有1個梯級流域,其中包含7個干流水電站及2個支流水電站共23臺水電機組,該梯級水電站結構如圖4所示。年負荷曲線取自該系統真實數據,如圖5所示。水火電機組的檢修時長均設置為實際檢修時長,如表4所示。T表示火電機組,H表示水電機組。

表4 水火電機組檢修時長(天)Tab.4 Maintenance time length of hydrothermal units (day)

圖4 水電站之間的梯級耦合關系Fig.4 Cascade coupling relationship between hydropower stations

圖5 系統年負荷曲線Fig.5 Annual load curve

算例設計為使用SAVLR方法求解與使用CPLEX求解器直接求解時求解結果的對比,以驗證SAVLR方法在大規模系統上的適用性與有效性。計算結果如表5。

表5 省級實際水火電1151節點86機系統中采用SAVLR和CPLEX獲得的結果對比Tab.5 Comparison of results obtained through SAVLR and CPLEX in real provincial 1151-bus and 86-unit system

由表5可以發現,采用CPLEX求解器直接求解花費的時間為3h34min4s,而使用SAVLR法求解時,若引入的對耦合約束的懲罰項系數γ取20,則計算在4次迭代花費34min56s后收斂,求解結果相對直接求解得到的結果誤差約為1.65%,在該次計算中,共求解火電子問題5次,水電子問題3次。而γ取50時,在3次迭代花費23min35s后收斂,相對誤差約為4.13%,在該次計算中,共求解火電子問題3次,水電子問題2次。但采用ADMM求解時,迭代過程不收斂。

從求解結果和計算時間可以看到,SAVLR法極大地減少了求解問題所需要的時間,求解得到的結果與直接求解得到的結果偏差不大,證明了SAVLR方法在大規模檢修問題上的有效性。同時,子問題計算數量的不同也說明了算法的一個優勢,在每次迭代中無須求解所有的子問題以更新拉格朗日乘子,這一點大大減少了計算時間。并且隨著γ取值的增大,計算時間和迭代次數會減少,但結果的相對誤差同樣會有所增大。當然,若γ取得太大則會造成過早收斂,最終結果誤差過大的問題,是不可取的。

除了算法的有效性以外,模型的合理性同樣需要說明,圖6和圖7分別為目標函數中w取0時優化得出的水電機組和火電機組的檢修計劃安排,圖中深色塊表示檢修,淺色塊表示不檢修,橫坐標為日期,縱坐標為機組編號。

圖6 火電機組的檢修計劃Fig.6 Maintenance schedule of thermal units

如圖7所示,由于目標函數中考慮了棄水量,優化得出水電機組的檢修計劃基本上分布在自然來水較為豐富的汛期之外。與此同時,火電機組因為功率平衡及備用約束,檢修計劃大部分安排在汛期附近,也就是說,火電機組為照顧水電機組在汛期盡量多發滿發而檢修,水火電機組間的檢修計劃大致上呈現錯峰的特點。

圖7 水電機組的檢修計劃Fig.7 Maintenance schedule of hydro units

以上計算結果均為目標函數中的懲罰系數w設置為0計算得出。以下設置懲罰系數分別為106和206分別進行計算,計算結果如表6所示。

表6 不同懲罰系數下計算結果對比Tab.6 Comparison of results corresponding to different penalty coefficients

表6中,懲罰系數w=0時的計算結果顯示,所有機組的檢修計劃均被調整。而懲罰系數w=106時,檢修計劃被調整的機組數由w=0時的71臺減少到35臺,但相應的火電運行成本上升了0.437%,棄水成本上升了5.85%;當w=206時,調整機組數進一步減少到27臺,而火電運行成本和棄水成本相較w=0時分別上升了0.912%和8.47%。顯然,當懲罰系數w增大時,被調整檢修計劃的機組數減少,但火電運行成本和棄水成本也會相應地增加。在實際中,電力調度中心可以根據對發電廠滿意度的重視程度靈活調整該系數,若設置為0,則求解得出經濟性最優的檢修計劃;若設為較大的值,則以犧牲一些經濟性為代價在一定程度上照顧發電廠的檢修意愿。

4 結論

本文以天為時間單位,建立了水火電機組年度檢修計劃MILP模型,綜合考慮了水火電機組的運行要求及發電廠的滿意度,并且包含了梯級水電耦合及線路傳輸容量限制等復雜約束。為加快問題求解速度,提出將原問題分解為水電和火電兩個子問題,并采用SAVLR法進行求解,該算法具有收斂性能好、計算速度快的優點。

本文基于3節點6機系統算例說明了求解方法的有效性,通過實際1151節點86機系統算例驗證了該方法的有效性,為求解大規模機組檢修計劃優化問題提供了新的求解思路。同時,通過實際系統對建立的模型進行了驗證,證明該模型能夠合理地安排水火電機組檢修任務的開始時間,使水電機組在汛期能夠盡量滿發;同時,電網公司可通過調整目標函數中第二項系數w的大小來衡量對發電廠檢修意愿的傾向程度,具有實際工程價值。

主站蜘蛛池模板: 久久成人免费| 亚洲成a人片| 欧美日韩精品一区二区视频| 一本一道波多野结衣av黑人在线| 婷婷色中文网| 99这里精品| 伊人激情久久综合中文字幕| 91黄色在线观看| 黄色成年视频| 国产凹凸一区在线观看视频| 91成人在线观看视频| 亚洲视频一区| 久久香蕉国产线看观看精品蕉| 国产成人午夜福利免费无码r| 欧美影院久久| 久久鸭综合久久国产| 中文字幕人妻无码系列第三区| 国产在线一二三区| 国国产a国产片免费麻豆| 国产又色又刺激高潮免费看| 成人日韩视频| 香蕉色综合| 毛片视频网| 日本免费高清一区| 91小视频版在线观看www| 国产95在线 | 91久久国产成人免费观看| 免费啪啪网址| 97视频免费看| 一级毛片基地| 伊人色婷婷| 亚洲中久无码永久在线观看软件| 国产办公室秘书无码精品| 久久久国产精品免费视频| 欧美日韩一区二区三| 波多野结衣一区二区三区四区| 久青草免费在线视频| 亚洲天堂成人在线观看| 天天综合色天天综合网| 亚洲六月丁香六月婷婷蜜芽| 69视频国产| 欧美性猛交一区二区三区| 强乱中文字幕在线播放不卡| 午夜国产在线观看| 久久99国产综合精品1| aaa国产一级毛片| 狂欢视频在线观看不卡| 亚洲天堂日本| 国产综合在线观看视频| 亚欧乱色视频网站大全| 露脸真实国语乱在线观看| 亚洲精品福利网站| 无码视频国产精品一区二区| 亚洲AV无码乱码在线观看代蜜桃 | 天堂成人在线| 2020精品极品国产色在线观看 | 亚洲中久无码永久在线观看软件| 中文无码精品a∨在线观看| 亚洲中字无码AV电影在线观看| 久久久久九九精品影院| 亚洲精品va| 一级爆乳无码av| 夜色爽爽影院18禁妓女影院| 啊嗯不日本网站| 国产成人91精品| www欧美在线观看| 国产黄色免费看| 午夜精品区| 波多野结衣在线se| 永久免费精品视频| 亚洲天堂免费观看| 伊人久久大线影院首页| 青青草欧美| 最新痴汉在线无码AV| 99久久无色码中文字幕| 五月婷婷丁香综合| 精品久久久久成人码免费动漫 | 亚洲专区一区二区在线观看| 欧美国产另类| 精品无码视频在线观看| 99久久99这里只有免费的精品| 19国产精品麻豆免费观看|