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

電-碳-綠證市場耦合下發電商報價與出清雙層優化

2024-03-06 09:18:00陳荃張丹宏鄭淇源郇嘉嘉趙敏彤朱建全
南方電網技術 2024年1期
關鍵詞:優化模型

陳荃,張丹宏,鄭淇源,郇嘉嘉,趙敏彤,朱建全

(1.廣東電網有限責任公司,廣州 510060;2.華南理工大學電力學院,廣州 510640)

0 引言

隨著“碳達峰、碳中和”戰略目標的提出[1],我國先后啟動了碳交易市場和綠證交易市場[2],促進碳減排和可再生能源的發展。與此同時,碳交易市場和綠證交易市場的引入也對電力市場產生了深刻影響。由于碳成本被納入發電商報價決策中,碳市場將直接影響電力市場的出清順序,改變現有發電側利益格局[3]。另一方面,綠證市場為可再生能源發展提供了額外的資金支持。可再生能源發電商通過出售綠證獲取收益,提高了市場競爭力。碳市場、綠證市場與電力市場間存在著復雜的相互作用關系,對發電商的經營決策和利益產生了重大影響。在此背景下,分析市場主體策略行為及電力市場均衡有助于發電商準確判斷市場運行狀態,從而獲取自身最優報價策略。

目前已有較多關于電力市場、綠證市場和碳市場協同作用及市場主體策略行為的研究。在模型方面,已有的研究可以分為系統動力學模型和一般均衡模型兩類。文獻[4-6]利用系統動力學方法建立了碳市場、綠證市場和電能量市場交互作用分析模型。這種模型相對簡單,但由于忽略了市場中各主體追求利益最大化的行為,往往難以幫助市場參與者分析獲取收益最大化的策略。同時,該模型依賴于大量參數設置,在參數設置不合理時,難以準確模擬真實市場的運行。為了解決這個問題,部分研究采用一般均衡模型對不同市場主體的報價策略進行建模。文獻[7]基于經典的古諾博弈競爭理論構建了電力市場的一般均衡模型。但由于古諾模型需要對電力市場進行大幅簡化(如沒有考慮輸電網絡等),所得結果可能與實際電力市場出現較大偏差。文獻[8]建立了計及潮流約束的雙層優化模型。其中,上層模型為傳統發電商和可再生能源發電商的策略報價;下層模型為計及綠證市場作用的日前電能量市場與實時電能量市場出清。但是,該文獻忽略了碳減排等因素,難以適用于電力市場和碳市場的協同分析。為此,文獻[9]建立了多時間耦合的電-碳市場分析模型,模擬了電力市場和碳市場之間的交互影響,能夠幫助傳統發電商在電-碳耦合市場中進一步優化決策。但是,該模型未考慮碳市場對不同時間尺度的電能交易的影響。

在算法方面,基于一般均衡理論的發電商決策模型是一個多時段的非線性雙層規劃問題,求解較為困難。文獻[10-11]基于罰函數方法將雙層優化問題轉換為無約束優化問題進行求解。但是,該方法依賴于懲罰因子等參數的設置,當參數設置不合理時,算法的收斂速率較慢。文獻[12-14]基于KKT(Karush-Kuhn-Tucker)條件將雙層優化模型轉換為單層優化模型,并通過強對偶定理和二進制拓展法進一步將單層優化模型轉換為混合整數線性規劃,最后調用求解器Gurobi 或CPLEX 進行求解。但是,由于轉換后的模型含有大量的整數變量,帶來的計算負擔較大。為了解決這個問題,文獻[15]提出了一種基于進化理論的優化算法,通過逼近雙層模型中的下層模型值函數減少了計算時間。文獻[16]提出一種基于嵌套理論和Kriging 插值的求解算法。該算法通過嵌套進化方法產生樣本,并用Kriging 插值近似下層模型的最優決策,能夠有效求解含有非線性約束的雙層模型。然而,這種近似化的處理方式所得結果與實際情況存在一定的偏差。

在這種背景下,本文首先介紹了電力、綠證和碳市場的交易框架,并在此基礎上建立了電-碳-綠證耦合的發電商策略報價與出清模型。其次,將該模型轉換為單層優化問題,并提出值函數近似算法對其進行求解。與現有方法相比,本文的主要創新點包括:1) 現有的電-碳-證耦合的交易模型一般采用階梯報價方式,而本文采用的是分段線性報價方式,能更好地反映發電商的邊際成本,并支撐發電商更為靈活地調整其報價策略;2) 現有的雙層優化方法主要基于KKT 條件等將雙層優化模型轉換為單層優化模型,而本文首次建立一個最優值函數以反映上下層決策之間的關系,不需要引入新的整數變量便可將雙層優化模型轉換為單層優化模型,從而實現電-碳-綠證耦合下報價與出清雙層模型的快速、有效求解。

1 市場架構

1.1 市場交易框架

圖1 給出了引入綠證交易、碳交易后的市場交易框架。其中,可再生能源發電商與傳統發電商共同參與日前、實時電能量市場的投標競爭,以獲得中標出力。由于風電、光伏發電等已納入綠證的核發對象,可再生能源發電商可通過參與綠證市場獲得額外的收益。同時,傳統發電商通過參與碳市場購買或出售碳配額差額,完成履約期內的碳配額繳納。電力用戶則通過參與日前、實時電能量市場購買電力。

圖1 市場交易框架Fig.1 Trading framework of markets

1.2 市場模型假設

1.2.1 電能量市場

電能量市場包括日前市場和實時市場。交易主體包含可再生能源發電商、傳統發電商和電力用戶。在日前和實時市場中發電商需提交價格-容量關系曲線。本模型中發電商采用分段線性方式報價,且各分段斜率單調遞增[17],如圖2所示。

圖2 分段線性報價曲線Fig.2 Piecewise linear bidding curves

考慮到大多數發電商的運行成本為二次函數,邊際成本則為一次函數[18]。因此,線性報價方式有利于發電商以邊際成本函數為基礎進行報價,符合各市場主體的利益。同時,在該報價方式下量價曲線的截距和各分段斜率均可調整,發電商的報價更為靈活,有利于電力市場的充分競爭。

日前和實時市場出清的線路模型為考慮傳輸功率限制的直流潮流模型[13]。同時,對于可再生能源機組出力的不確定性與波動性,本文參照文獻[19-20]設定可再生能源機組的出力預測誤差滿足標準高斯分布,并通過蒙特卡洛抽樣獲得實時市場的多個場景,以模擬發電機組實際出力與預測出力的偏差[8]。

1.2.2 綠證市場

綠證市場的交易主體為可再生能源發電商、電力零售商等。可再生能源發電商通過出售綠證獲得額外的收益。電力零售商等通過購買綠證完成配額指標。本文參照文獻[8]設定了綠證的最高價格為配額義務主體無法完成配額任務時的罰款價格。同時,綠證交易的量價關系滿足古諾產量競爭模型[8]。由于可再生能源發電商出售的綠證數量在電能量市場出清后方可確定,本文的模型假定綠證市場的交易發生在日前和實時市場之后,且按照獨立的交易時間尺度開展。

1.2.3 碳市場

碳市場的交易主體為傳統發電商。傳統發電商需在履約期結束時繳納該年度碳排放量對應的碳配額。在本模型中假定傳統發電商將初始免費發放的碳配額分解到較短的時間尺度,并結合實際發電量對應的碳排放量最終確定在碳市場中的交易量。將配額分解納入碳市場模型中可以更精準地刻畫傳統發電商的策略行為,從而更好地模擬碳市場的運行過程[9]。在交易時序上,本文假定碳市場的交易在日前和實時市場之后開展,交易的時間尺度與綠證市場一樣。同時,模型中的碳價由外部輸入,不考慮碳市場的具體出清過程。

1.2.4 各市場銜接關系

在碳市場中,傳統發電商通過碳交易優化配置碳資源以成本效益最優的方式完成碳配額考核。在綠證市場中,可再生能源發電商通過綠證交易獲取綠色電力的環境價值收益。同時,兩類發電主體共同參與電力市場競爭,獲取電能量價值對應的收益。綠證市場和碳市場也因發電主體在電力市場中的競爭關系而相互耦合在一起。需要說明的是,目前初始碳配額免費發放、可再生能源消納責任權重的考核要求較低,導致各市場間的耦合作用相對較弱。未來,隨著碳配額考核和可再生能源消納責任權重考核的進一步收緊,各市場間的關聯性將進一步增強。

2 多主體優化決策與市場出清模型

本文所考慮的多主體優化決策及市場出清模型是一個雙層博弈問題,其結構如圖3 所示。在上層模型中,各主體通過最大化其利潤確定自身的報價曲線,并將其傳遞給下層;在下層模型中,日前電能量市場、實時電能量市場、綠證市場分別出清,得到節點邊際電價和各主體的出清電量后將其返回給上層模型。

圖3 雙層多主體優化決策模型Fig.3 Bi-level optimization decision model of multi-entity

2.1 上層模型

2.1.1 目標函數

以發電商j在各市場中的總利潤最大化為目標,有:

2.1.2 約束條件

1) 報價約束

在本文中,發電商采用分段線性方式報價,且各分段斜率單調遞增,滿足如下約束。

2) 綠證出售數量約束

可再生能源發電商j在場景ω下出售的綠證數量滿足以下約束。

2.2 下層模型

2.2.1 下層模型1:日前市場出清

2.2.1.1 目標函數

以最小化日前市場的負社會福利(即最大化日前市場的社會福利)為目標,有:

式中:d為負荷的編號;k為負荷的報價容量段的序號;γi,b,t為發電機組i在時段t的報價曲線中第b分段的起點對應的縱坐標值(價格);λd,k,t為負荷d在時段t、容量段k的效用值;為日前市場中負荷d在時段t、容量段k的中標容量;為發電機組i在容量段b的出力上限。同時,社會總福利值由消費者盈余減去生產者盈余得到。

2.2.1.2 約束條件

日前市場出清需滿足電網、機組、負荷等約束,具體如下。

1) 直流潮流約束

2) 節點相角約束

2.2.2 下層模型2:實時市場出清

2.2.2.1 目標函數

以最小化實時市場的負社會福利(即最大化實時市場的社會福利)為目標,有:

2.2.2.2 約束條件

實時市場出清也需滿足電網、機組、負荷等約束,具體如下。

1) 直流潮流約束

5) 負荷功率約束為確保負荷各分段功率不超過上下限,需滿足以下約束。

2.2.3 下層模型3:綠證市場的出清

綠證市場的出清模型可通過式(29)—(32)描述。

式(29)建立了古諾模型中的綠證逆需求線性函數;約束(30)規定綠證的最高價格為配額義務主體無法完成配額任務時的罰款價格;約束(31)為計算逆需求線性函數的斜率,可由配額比例和歷史數據確定;約束(32)用于計算電力用戶的最大綠證需求量。

3 值函數近似算法

從數學上看,上述的模型是一個多時段的非線性雙層規劃問題。傳統算法主要通過KKT 條件將雙層問題轉換為單層問題,進而采用二進制拓展法處理單層問題中的互補松弛約束,問題求解困難且耗時較長。因此,本文提出一種基于值函數近似的雙層模型求解算法以降低計算復雜度,提高計算效率。

3.1 基本思路

首先,建立上層決策變量與下層目標函數的映射關系,即最優值函數。其次,基于多項式基建立函數的近似結構,并用最小二乘法求解待定系數。然后,基于近似值函數將雙層模型轉換為單層模型。最后,采用對角化算法求解個體間的相互博弈問題。

3.2 構建最優值函數

一般的雙層優化問題可表示為:

式中:x為上層優化問題的決策變量;y為下層優化問題的決策變量;p為下層優化問題中約束的序號;P為下層優化問題中約束的數量;k為上層優化問題中約束的序號;K為上層優化問題中約束的數量;Gk(x,y)為 上 層 優 化 問 題 中 第k條 約 束;gp(x,y)為下層優化問題中第p條約束;F(x,y)為上層優化問題的目標函數;f(x,y)為下層優化問題的目標函數。

為了替代下層模型,引入映射函數來表示上層優化問題的決策變量與下層優化問題的最優值函數之間的映射關系,該關系可以用圖4 來表示。其數學表達式如式(36)所示。

圖4 值函數映射關系Fig.4 Mapping relationship of the value function

式中:V為最優值函數;φ(?)為映射函數。

在本模型中,下層模型最優值函數的輸入為各發電商的報價曲線,輸出為市場的總社會福利值或節點邊際電價,具體可表示為式(37)。

式中αR為最優值函數的輸入變量(即發電商報價曲線的參數信息)。

3.3 基于多項式基的值函數近似技術

根據文獻[21-22],值函數可表示為關于輸入變量αR的多項式基函數的線性組合,即:

式中:φ?(αR)為值函數V的近似值;l為多項式基函數的序號;dl為第l個多項式基函數的系數;?l(αR)為第l個多項式基函數;F 為所有多項式基函數構成的集合,具體可表示為式(39)。

選擇更高一階多項式基函數的零點作為采樣點,則采樣點的數量為:

式中N為采樣點的總數量。

將采樣點及其對應的值函數觀測值代入式(38),可得到關于系數dl的回歸方程為:

式中:?N,Ns為第Ns個多項式基函數在第N個采樣點處的值;dNs為第Ns個多項式基函數的系數;φN為φ(αR)在第N個采樣點處的值函數值,可通過調用求解器解下層模型得到。

根據文獻[23-24],利用最小二乘法求解多項式基函數的待定系數dl,有:

式中:D為待定系數dl構成的矩陣;Φ為多項式基函數?N,Ns構成的矩陣;ψ為值函數觀測值φN構成的矩陣。

3.4 雙層模型轉換為單層模型

對于發電商j,對應的單層模型可表示為:

目標函數式(43)為最小化發電商j的負利潤(最大化發電商j的利潤);約束式(44)用于確保日前電能量市場的負社會福利最小;約束式(45)用于計算日前電能量市場在時段t的節點邊際電價;約束式(46)用于確保實時電能量市場的負社會福利最小;約束式(47)用于計算實時電能量市場在場景ω、時段t的節點邊際電價。

此外,在單層模型的求解過程中還需滿足約束式(2)—(11)、約束式(13)—(20)和約束式(22)—(32)。

3.5 基于對角化算法的多個體博弈求解

將雙層模型轉換為單層模型后,各發電商的策略報價及出清問題可描述為多主體博弈問題。本文采用文獻[13]提出的對角化算法對其進行求解。對角化算法是通過迭代求解的方式得到各主體的最優報價曲線。當求解某一個體的最優報價曲線時其他個體的報價曲線采用上一輪迭代求解的結果。同時,電力用戶的報價在迭代過程中始終保持不變。算法的終止條件是所有個體本輪報價與前一輪報價的差值小于某一給定值。將對角化算法用于求解本文所提的多主體博弈問題的具體流程如圖5所示。

圖5 對角化算法流程圖Fig.5 Flow chart of diagonalization algorithm

4 算例分析

4.1 基礎數據

為驗證所提算法的有效性,采用一個經修改的IEEE 30節點的系統進行仿真分析。該系統包含6臺發電機組和11個負荷。其中,G1、G2、G5、G6分屬于傳統發電商1—4,相關參數依據文獻[25];G3、G4分屬于可再生能源發電商1和2,相關參數依據文獻[26-27]。同時,參照文獻[8]和[12],設置4個典型場景和6個仿真時段。其中,各場景及各時段的負荷和可再生能源機組的出力如圖6—9所示。

圖6 場景1下負荷和可再生能源機組出力Fig.6 Loads and outputs of renewable energy units in scenario 1

圖7 場景2下負荷和可再生能源機組出力Fig.7 Loads and outputs of renewable energy units in scenario 2

圖8 場景3下負荷和可再生能源機組出力Fig.8 Loads and outputs of renewable energy units in scenario 3

圖9 場景4下負荷和可再生能源機組出力Fig.9 Loads and outputs of renewable energy units in scenario 4

綠證市場中配額比例ρTGC0為30%,價格常數θTGC0為0.3,罰款價格λTGC0為200 元/MWh。碳市場中,碳交易價格λCT為40 元/t,碳排放因子EGi為0.7。

本文所有的仿真分析均在MATLAB 和GAMS平臺編程實現。計算機配置為:Intel Core-i7 處理器,主頻為2.1 GHz,內存為16 GB。

4.2 算法性能測試

4.2.1 值函數近似精度測試

本文采用多項式基函數和最小二乘法對下層模型的最優值函數進行近似,進而將電力市場雙層模型轉換為單層模型。為了驗證所提算法的準確性,以下將比較近似前后的下層模型最優值函數。

首先,利用蒙特卡洛方法采樣生成100 個測試點,每個測試點均包含發電商的報價信息。其次,調用GAMS 中CONOPT 求解器求解下層模型,獲得最優值函數的真實值。最后,根據式(38)計算最優值函數的近似值。

圖10—11 分別為不同測試點下利用原始下層模型及其近似值函數輸出的社會福利值。其中,前后者分別被視為真實值和近似值以評估值函數的準確性。從圖10—11 可見,最優值函數的近似值與真實值相對誤差小于1%。這表明所提的值函數近似算法具有較高的準確性,能夠滿足雙層模型求解精度的要求。

圖10 最優值函數的近似值與真實值Fig.10 Approximate values and true values of the optimal value function

圖11 最優值函數的近似誤差Fig.11 Approximation errors of the optimal value function

4.2.2 算法效率測試

為了測試所提算法的計算效率,將報價允許誤差(收斂條件)ε從0.02 增加至0.10,以步長0.02遞增,其他參數保持不變。不同收斂條件下所提算法的運行時間如表1所示。

表1 不同收斂條件下算法運行時間Tab.1 Runtimes of the algorithm under different convergence conditions

從表1 可見,隨著報價允許誤差的增大算法的計算時間呈現下降趨勢,且最大運行時間為128 s,平均運行時間為77 s。因此,本文所提算法具有較高的求解效率,能夠滿足分析實際電力市場的需求。

4.3 運行結果分析

4.3.1 各發電商的中標出力分析

按本文所提方法對各發電商的報價策略進行優化,得到各發電商的報價曲線和各發電機組在市場中的中標出力,如圖12—13 所示。考慮到目前現貨市場普遍采用階梯形式進行申報,在優化得到線性報價曲線后,將線性報價曲線還原為階梯報價曲線。

圖12 時段1各發電商的報價曲線Fig.12 Bidding curves of each power producer during period 1

圖13 時段1—6各機組的中標出力Fig.13 Bid-winning outputs of each unit during period 1—6

從圖12—13 可見,傳統發電商3(對應機組G5)的報價最高,因而在實時市場中獲得的中標出力最少。雖然可再生能源發電商2(對應機組G4)在容量段1(0—100 MW)的報價最低,但由于可再生能源機組出力的不確定性,發電機組在實時市場中獲得的中標出力也較小。相對而言,傳統發電商1(對應機組G1)和傳統發電商4(對應機組G6)的報價相對較低,在實時市場中獲得較多的中標出力,具備較強的市場競爭力。

4.3.2 各典型場景下的出清電價分析

在前文所述的各種典型場景下利用所提方法制定各發電商的報價策略。當市場達到均衡時根據各發電商的報價曲線及機組中標出力計算實時市場系統邊際電價,結果如圖14所示。

圖14 場景1—4下各時段實時出清電價Fig.14 Real-time clearing prices in each period of scenarios 1—4

從圖14 可見,電力的供需平衡關系對市場出清電價具有較大影響。其中,時段4 的市場出清電價最高,為143~163元/MWh。這是因為時段4處于用電高峰期,負荷達到了870 MW。時段1和時段2的負荷需求較小,市場出清電價也明顯降低,僅為101~107元/MWh。另一方面,場景2下的市場出清電價明顯高于其他場景下市場出清電價,這是因為場景2 下可再生能源機組的出力最低。這種供電緊張的局面造成了市場出清電價顯著上升。

4.3.3 碳交易價格對市場均衡的影響

為了探究碳交易價格λCT對電力市場均衡的影響,將碳交易價格λCT從40 元/t 提高至80 元/t,以步長10元/t遞增,并分別記為碳價1—碳價5,其他參數與4.3.1 小節相同。運行所提的發電商策略報價與出清方法,結果如圖15所示。

圖15 碳交易價格對實時出清電價的影響Fig.15 Influence of carbon trading prices on real-time clearing price

從圖15 可見,隨著碳交易價格的上升實時市場的出清電價呈現上升趨勢。這是因為傳統發電商的碳排放成本增大,其提交的報價曲線將相應上移,市場出清電價因而提高。其中,時段3 和時段4 的負荷最大,傳統機組的出力最多,因而市場出清電價受碳交易價格的影響也最大,呈現顯著上升趨勢。

4.3.4 可再生能源滲透率對市場均衡的影響

為了探究可再生能源滲透率對市場均衡的影響,將可再生能源滲透率從25%提高至45%,以步長10%遞增,其他參數與4.3.1 小節相同。運行所提的發電商策略報價與出清方法,市場運行結果如圖16所示。

圖16 可再生能源滲透率對實時出清電價的影響Fig.16 Influence of renewable energy penetration rates on realtime clearing prices

從圖16 可見,隨著可再生能源滲透率的增長,實時市場的出清電價呈現下降趨勢。這是因為可再生能源發電商運行成本較低,且可通過銷售綠證獲取額外收益,其在電能量市場上的報價曲線往往低于傳統發電商。因此,可再生能源占據的市場份額越大市場的出清電價就越低。

4.3.5 不同市場耦合的效果分析

為了探究不同耦合市場對可再生能源消納和社會經濟效益的影響,分別在電力市場、電力-碳耦合市場、電力-綠證耦合市場及電力-碳-綠證耦合市場下進行仿真分析,結果如表2所示。

表2 不同市場下可再生能源發電占比及經濟效益Tab.2 Proportions of renewable energy generation and economic benefits under different markets

從經濟性的角度看,電力-綠證耦合市場的社會經濟效益最高,電力-碳耦合市場的社會經濟效益最低。這是因為,引入綠證市場后可再生能源發電商可通過銷售綠證獲取額外收益,其在電力市場的報價降低,社會經濟效益也就相應增大。而引入碳市場后傳統發電商的發電成本增大,致使其在電力市場中提高報價,造成社會經濟效益減小。考慮可再生發電商在綠證市場的銷售收入及傳統發電商在碳市場額外支出的成本后電力-碳-綠證耦合市場的社會經濟效益介于電力-綠證耦合市場和電力-碳耦合市場之間。

從可再生能源消納的角度看,電力-碳-綠證耦合市場下可再生能源發電占比最高。這是因為,引入綠證市場與碳市場后可再生能源發電商的報價降低,而傳統發電商的報價提升,這提高了可再生能源機組在市場出清時的中標出力。因此,電力-碳-綠證耦合市場能夠促進可再生能源發展,并兼顧經濟效益,總體效果較好。

5 結論

本文首先建立了電能量市場、綠證市場和碳市場耦合的多主體報價與出清雙層優化模型,并提出了一種基于值函數近似的優化算法對其求解,得到主要結論如下。

1) 所提算法具有較高的求解精度和計算效率。在精度方面,近似值函數的輸出值與真實值幾乎一致,能夠準確描述下層模型中發電商的報價對社會福利及邊際電價的響應特性。在速度方面,傳統方法難以求解本文模型,而所提方法在較短時間內即可獲得近似最優解。

2) 負荷需求和可再生能源發電對市場出清電價具有較大的影響。當負荷需求較大或可再生能源機組出力較低時市場的出清電價較高;隨著可再生能源滲透率的增長,市場的出清電價呈現下降趨勢。

3) 隨著碳交易價格的上升,市場的出清電價呈現上升趨勢,且傳統機組的出力越多,市場出清電價受碳交易價格的影響越明顯。

4) 電力-碳-綠證耦合市場能夠促進可再生能源發展,并兼顧經濟性,總體社會效益較好。

在后續研究中將進一步考慮發電機組開停及其成本對電-碳-證市場耦合下發電商報價與市場出清的影響,以更好地促進新能源消納與電網安全穩定運行。

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产欧美日本在线观看| 免费a级毛片18以上观看精品| 国产视频a| h网址在线观看| 亚洲精品视频网| 四虎成人精品| 91小视频版在线观看www| 91小视频在线观看免费版高清| 亚洲人成色在线观看| 美女视频黄又黄又免费高清| 一级毛片免费不卡在线视频| 久操中文在线| www.99在线观看| 奇米影视狠狠精品7777| 激情综合图区| 免费国产高清精品一区在线| 91精品国产情侣高潮露脸| 91在线日韩在线播放| 亚洲无码高清视频在线观看| 亚洲精品无码久久久久苍井空| 四虎成人免费毛片| 香蕉久人久人青草青草| 国产丝袜第一页| 亚洲无码精品在线播放| 国产精品亚欧美一区二区| 久草视频精品| 色哟哟精品无码网站在线播放视频| 国产精选自拍| 久草视频一区| 久久精品国产亚洲AV忘忧草18| 久久99这里精品8国产| 国产免费怡红院视频| 亚洲国产91人成在线| 国产精品成人免费视频99| 朝桐光一区二区| 97免费在线观看视频| 亚洲成人手机在线| 成人国产免费| 日韩欧美中文| 91在线丝袜| 超碰91免费人妻| 国产欧美中文字幕| 日韩在线第三页| 日韩精品中文字幕一区三区| 欧美另类第一页| 欧美福利在线| 99久久这里只精品麻豆| 国产精品任我爽爆在线播放6080| 黄色片中文字幕| 欧美自慰一级看片免费| 亚洲无码高清免费视频亚洲| 99re这里只有国产中文精品国产精品 | 精品黑人一区二区三区| 老熟妇喷水一区二区三区| 97色婷婷成人综合在线观看| 女人18毛片久久| 国产主播在线观看| 中国黄色一级视频| 免费高清自慰一区二区三区| 天堂网亚洲系列亚洲系列| 91国内视频在线观看| 欧美性色综合网| 欧美一区福利| 国产69囗曝护士吞精在线视频| 国产一区二区人大臿蕉香蕉| 亚洲国语自产一区第二页| 丰满人妻一区二区三区视频| 亚洲二区视频| 538国产视频| 亚洲国产看片基地久久1024| 91人人妻人人做人人爽男同| 亚洲第一成人在线| 青青青国产免费线在| 91毛片网| 日韩精品一区二区三区中文无码| 国产成人麻豆精品| 91福利在线看| 久热中文字幕在线观看| 免费全部高H视频无码无遮掩| 亚洲国产欧洲精品路线久久| www精品久久| 国产美女无遮挡免费视频网站|