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

基于二階錐規(guī)劃凸松弛的三相交直流混合主動配電網(wǎng)最優(yōu)潮流

2021-05-12 06:24:48巨云濤張若思
電工技術(shù)學(xué)報 2021年9期
關(guān)鍵詞:模型

巨云濤 黃 炎 張若思

(中國農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院 北京 100083)

0 引言

近幾年,隨著分布式電源廣泛接入配電網(wǎng),傳統(tǒng)的被動配電網(wǎng)轉(zhuǎn)變?yōu)檎{(diào)度資源豐富的主動配電網(wǎng)[1-2]。

20 世紀(jì)60 年代提出的最優(yōu)潮流(Optimal Power Flow, OPF)是基于數(shù)學(xué)的最優(yōu)化理論[3]。OPF 是一個非凸非線性問題,其求解思路一般是將非凸問題進(jìn)行凸松弛轉(zhuǎn)換成凸問題以保證可靠求解[4-5]。基于凸優(yōu)化理論的凸規(guī)劃不但保證了計算結(jié)果為全局最優(yōu)解,還能夠為非凸問題提供較好的上下邊界[6-7]。

在電力系統(tǒng)中應(yīng)用最為廣泛的凸優(yōu)化方法是二階錐規(guī)劃(Second Order Cone Programming, SOCP)和半正定規(guī)劃(Semi-Definite Programming, SDP),以及近幾年提出的二次凸松弛技術(shù)(Quadratic Convex Relaxation, QC Relaxation),前兩者的求解方法大多是基于內(nèi)點法[8]。文獻(xiàn)[9-10]中首次提出了SOCP 模型在環(huán)網(wǎng)和輻射網(wǎng)潮流問題中的應(yīng)用,并且驗證了模型的有效性,模型結(jié)果與牛頓法的求解結(jié)果非常相近。文獻(xiàn)[11]建立了SOCP 的拓展模型,以解決環(huán)網(wǎng)中電壓的相位差無法應(yīng)用SOCP 約束的問題,這一變化是錐理論從潮流問題跨入最優(yōu)潮流問題的重要一步。文獻(xiàn)[12]提出了考慮電壓穩(wěn)定約束的最優(yōu)潮流模型,將SOCP 松弛應(yīng)用于電壓穩(wěn)定控制中。文獻(xiàn)[13]構(gòu)建了以最優(yōu)能量流為目標(biāo)函數(shù)的電氣-天然氣綜合能源系統(tǒng)詳細(xì)模型,并采用二階錐規(guī)劃將非凸規(guī)劃問題轉(zhuǎn)換為凸規(guī)化問題求解。文獻(xiàn)[14-16]使用SOCP 松弛,系統(tǒng)性地將原支路最優(yōu)潮流模型轉(zhuǎn)換為SOCP,保證模型所獲得的解為全局最優(yōu)解。文獻(xiàn)[17-18]通過計算不同系統(tǒng)發(fā)現(xiàn)SOCP 對于大規(guī)模的電力系統(tǒng)有較高的求解效率。文獻(xiàn)[19]提出了一種基于SOCP 松弛的配電網(wǎng)動態(tài)最優(yōu)潮流模型,并討論了在接入不同配電網(wǎng)元件的情況下凸松弛模型的計算準(zhǔn)確性。文獻(xiàn)[20]建立了交直流混合系統(tǒng)最優(yōu)潮流的SOCP 模型,提出一種考慮可再生能源不確定性的可調(diào)魯棒優(yōu)化方法。文獻(xiàn)[21]建立了直流系統(tǒng)中的最優(yōu)潮流SOCP 松弛模型,并提出一種采用交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)的分布式計算方法。文獻(xiàn)[22]建立了配電網(wǎng)的三相無功優(yōu)化模型,并且將混合整數(shù)法與SOCP 結(jié)合,得到的松弛模型既滿足準(zhǔn)確性要求又具有較高的計算效率。雖然SDP 松弛的精確程度較高,但計算較為復(fù)雜,增加了問題的求解難度。隨著問題規(guī)模的增大,SOCP 在迭代計算上具有更大的優(yōu)勢[23]。

綜合考慮國內(nèi)外現(xiàn)有研究,與主動配電網(wǎng)相關(guān)的凸松弛研究較少,且現(xiàn)有研究并沒有給出分布式電源、有載調(diào)壓器等元件的詳細(xì)凸松弛模型。此外,現(xiàn)有關(guān)于Distflow 模型凸松弛的研究文獻(xiàn),并沒有對元件的物理特性進(jìn)行分析[24-25]。輸電網(wǎng)的凸松弛最優(yōu)潮流模型無法很好地應(yīng)用于主動配電網(wǎng),因此尋求更加準(zhǔn)確高效的主動配電網(wǎng)最優(yōu)潮流計算方法十分重要[26-27]。

本文首先提出含三相分布式電源、不同接線方式的有載調(diào)壓器、交直流混合系統(tǒng)換流器等模型的主動配電網(wǎng)SOCP 松弛模型,相比于Distflow 模型的凸松弛,本文考慮相間耦合,同時給出了分布式電源、有載調(diào)壓器等元件的物理特性;然后采用SOCP 松弛求解器求解,對比分析了不同求解器的計算效率差異。

1 模型介紹

由于一般最優(yōu)潮流模型約束中含有大量的非凸約束,不利于最優(yōu)模型求解,而引入SOCP 凸松弛方法,可將非凸約束轉(zhuǎn)換為凸約束條件,保證所得最優(yōu)解為全局最優(yōu)解。

1.1 普通三相支路

在一般潮流約束方程中,母線電壓、支路功率、支路電流一般滿足分別

式中,i、j分別為第i、j個三相母線;為母線i、j的三相電壓;Zij為母線i、j之間的三相線路阻抗;I˙ij為母線i流向母線j的電流。

在凸松弛的模型中,需要重新定義待求變量,重新定義的變量為二次型變量,與原模型的變量間存在一定的關(guān)系,但對于凸松弛模型,僅有重新定義的二次型變量為待求變量。定義二次型變量X為電壓、電流組成的矢量乘以該矢量的厄米爾特轉(zhuǎn)置所得到的厄米爾特矩陣為

式中,H 表示厄米爾特轉(zhuǎn)置,即先對矩陣取轉(zhuǎn)置,再對矩陣中的每一個元素取共軛;表示對復(fù)數(shù)相量取共軛;變量為厄爾米特矩陣,需滿足,其中Wi和lij為3 × 3矩陣。

對于重新定義的變量矩陣X,不僅需要滿足半正定約束,同時需要滿足秩1 約束,即變量矩陣的秩為1。秩1 約束是非凸約束,因此一般是將秩1約束松弛掉,即去掉秩1 約束進(jìn)行求解,并將所求結(jié)果代入秩1 約束進(jìn)行檢驗。若解滿足秩1 約束,則說明松弛是精確的,若解不滿足秩1 約束,則所求得的解無法還原出凸松弛前原模型的解,即所求解無意義。根據(jù)文獻(xiàn)[28-29],若變量矩陣X為厄米爾特矩陣且X≥0,則其所有的主子矩陣(principal submatrix)滿足行列式大于等于零的條件。因此將進(jìn)行SOCP 松弛后,需滿足矩陣X的所有2 ×2 的主子矩陣行列式大于等于零的條件,即

式(3)是一個2 ×2 厄爾米特矩陣,使用Sylvester準(zhǔn)則,式(3)可以轉(zhuǎn)換成,進(jìn)而可得到SOCP 的一般約束表達(dá)式

或SOCP 的標(biāo)準(zhǔn)約束形式

文獻(xiàn)[18]證明了上述轉(zhuǎn)換后的松弛條件與SOCP 松弛條件等價。

本文中的SOCP 約束應(yīng)轉(zhuǎn)換為

式中,(i,j)表示矩陣中的第i行,第j列。

經(jīng)過凸松弛后,式(1)變?yōu)?/p>

式(5)為式(2)中元素的線性組合。

1.2 分布式電源

分布式電源一般采用序分量模型,具有對稱的結(jié)構(gòu)參數(shù),因此線路中的分布式電源需滿足潮流約束。

其中

分布式電源接線圖如圖1 所示。

圖1 分布式電源接線Fig.1 Wiring diagram of distributed generation

這里,由于分布式電源采用三角形接線方式的調(diào)壓器隔離了零序電流,線路中零序電流為零,即

分布式電源如果不增加負(fù)序電流控制,負(fù)序電壓和電流關(guān)系符合阻抗特性,表示為阻抗方程,即

分布式電源的三相功率總和為常量,即

式中,Sg為分布式電源三相功率;表示對矩陣Sg的所有元素求和;St為已知三相功率總和。

式中,UG、IG、SG為厄爾米特矩陣變量,且。與1.1 節(jié)相似,松弛后的變量需滿足

式中,(i,j)為矩陣中的第i行,第j列。

定義“IG[1,1]表示矩陣IG中的第1 行,第1 列的元素”,則一般潮流約束(7)、(8)、(9)轉(zhuǎn)換為

1.3 直流部分

對于交直流混合系統(tǒng),也需要滿足基爾霍夫電流定律,以圖2 中的交直流混合系統(tǒng)為例,在交流區(qū)域1 中,對于母線2,有如下約束條件:

式中,1S交流區(qū)域1 中母線2 的三相注入功率;1PΔ為母線2 與直流支路節(jié)點1 之間的線路損耗;為計算直流支路損耗的系數(shù);Pc1為直流支路中節(jié)點1 的功率。

圖2 交直流混合系統(tǒng)電路拓?fù)鋱DFig.2 Circuit topology of AC-DC hybrid system

在直流區(qū)域中,直流支路滿足如下功率約束:

式中,Pc2為直流支路中節(jié)點2 的功率;Uc2為直流支路中節(jié)點2 的電壓;Ic為直流電流。

在交流區(qū)域2 中,對于母線3,滿足如下約束:

式中,2S為交流區(qū)域2 中母線3 的三相注入功率,2PΔ 為直流支路節(jié)點2 與母線3 之間的損耗;Pc2為直流支路中節(jié)點2 的功率。

在交流區(qū)域2 中必須將一個節(jié)點設(shè)定為參考節(jié)點來保證潮流系統(tǒng)的可解性,因此將母線3 設(shè)為恒定電壓源。

在直流部分,凸松弛的定義與普通三相支路相同,凸松弛后變量定義為

式中,Y為厄爾米特矩陣變量,且Y≥0;Wc1、Wc2為電壓二次型變量,且Wc1=Uc1Uc1,Wc2=Uc2Uc2;Sc1、Sc2為功率二次型變量,且Sc1=Uc1Ic,Sc2=Uc2Ic;lc為電流二次型變量,且lc=IcIc。變量Y應(yīng)滿是

式(15)~式(17)中的約束是變量矩陣(18)中元素的線性組合,對于部分一次項,需要做相應(yīng)的變換,如,電壓約束,可在等式兩邊同時乘上Ic,方程變換為該等式約束是變量矩陣(18)中元素的線性組合,且與原電壓約束等價。

1.4 有載調(diào)壓器

對于有載調(diào)壓器來說,其連接方式主要有Y 聯(lián)結(jié)和D 聯(lián)結(jié)兩種,這里以Y-Y 聯(lián)結(jié)和Y-D 聯(lián)結(jié)兩種聯(lián)結(jié)形式為代表。其中調(diào)壓器電壓比如果作為可控變量,會出現(xiàn)變量乘積的形式,即

式中,Um為調(diào)壓器一次電壓;Uj為調(diào)壓器二次電壓;y為調(diào)壓器電壓比。

需要注意的是,如果調(diào)壓器電壓比為未知,則為兩個變量相乘,由文獻(xiàn)[30]可知,式(20)可凸松弛為

式中,yl、yu和、分別為變量y和Uj的上、下界。

1.4.1 Y-Y 聯(lián)結(jié)

Y-Y 聯(lián)結(jié)如圖3 所示,在本模型中將調(diào)壓器的三相內(nèi)阻抗折算到一次側(cè),內(nèi)阻抗均滿足普通支路約束方程。

圖3 Y-Y 型聯(lián)結(jié)有載調(diào)壓器Fig.3 Y-Y connection group diagram of transformer

在一般潮流約束中,一、二次電壓、電流關(guān)系為

式中,S1、S2分別為一、二次側(cè)的支路功率。在一般潮流約束中,由電壓、電流約束可推出功率約束,但在凸松弛后的模型中,變量定義有所不同,電壓、電流約束并不能推出功率約束,故需增加有載調(diào)壓器一次側(cè)與二次側(cè)的功率約束方程,才能保證約束條件的完整性。

1.4.2 Y-D 聯(lián)結(jié)

Y-D 聯(lián)結(jié)如圖4 所示。在本模型中,將有載調(diào)壓器的三相內(nèi)阻抗折算到一次側(cè),內(nèi)阻抗均滿足普通支路約束方程。

圖4 Y-D 聯(lián)結(jié)有載調(diào)壓器Fig.4 Y-D connection group diagram of transformer

對于Y-D 型有載調(diào)壓器,在一般潮流約束中,一、二次電壓、電流關(guān)系為

其中

類似地,將式(24)進(jìn)行凸松弛后可轉(zhuǎn)換為

1.5 目標(biāo)函數(shù)

本文的目標(biāo)函數(shù)設(shè)為系統(tǒng)有功功率損耗,即網(wǎng)絡(luò)中有功功率損耗之和。

式中,ΔW為網(wǎng)絡(luò)有功功率損耗之和;式(26)中左邊項包括所有交流支路,即所有全相、兩相、單相運(yùn)行的交流支路;為3 × 3對角線矩陣,該矩陣對角線元素與列向量元素對應(yīng)行相等,非對角線元素值為零;Y ij為母線i和母線j之間的導(dǎo)納矩陣;式(26)右邊項包括所有直流線路,為第k條直流線路電流;為第k條直流線路電阻;Nc為直流線路總數(shù)。

2 算法分析

原始-對偶內(nèi)點法是求解SOCP 模型的有效方法之一,具有多項式時間復(fù)雜性,可以快速有效地求解大規(guī)模優(yōu)化問題[31]。原始-對偶內(nèi)點法對初始點的嚴(yán)格可行性要求較高,而實際問題中直接找到嚴(yán)格可行點并不容易,故改進(jìn)現(xiàn)有的原始-對偶內(nèi)點法對SOCP 的求解有著重要意義。光滑算法是近年來求解SOCP 問題的新方法,該方法是基于互補(bǔ)問題發(fā)展而來,不僅在理論上具有好的收斂性,而且實際計算效果也很好,也是求解SOCP 問題的常用方法之一[32]。

3 算例分析

3.1 IEEE 13 節(jié)點系統(tǒng)

為了驗證SOCP 模型的松弛精度,本文使用基于IEEE 13 節(jié)點的交直流混合系統(tǒng)進(jìn)行測試。系統(tǒng)結(jié)構(gòu)及參數(shù)詳細(xì)說明見文獻(xiàn)[33]。在測試算例中,將有載調(diào)壓器電壓比、分布式電源的功率等結(jié)果作為已知條件代入一般潮流算法中,所得到的一般潮流結(jié)果與凸松弛的優(yōu)化結(jié)果進(jìn)行比較。測試用的求解器是SCS。凸松弛模型的誤差分析主要有兩個指標(biāo):第一個是方均根誤差(Root Mean Squared Error,RMSE),其定義為預(yù)測值與真實值偏差的二次方和與觀測次數(shù)比值的二次方根,在本文中為二階錐規(guī)劃凸松弛所求各母線電壓結(jié)果與一般潮流各母線電壓結(jié)果的差的二次方和與母線個數(shù)的比值的二次方根;第二個是最大絕對誤差(Maximum Absolute Error, MAE),一般用來衡量絕對誤差的范圍,即二階錐規(guī)劃凸松弛所求各母線電壓結(jié)果與一般潮流各母線電壓結(jié)果的絕對誤差的最大值。兩個指標(biāo)相結(jié)合,合理表征凸松弛模型的計算精度。具體公式為

式中,N為系統(tǒng)節(jié)點總數(shù);為第i個母線所求電壓幅值;為第i個母線的電壓幅值。

3.1.1 交直流混合系統(tǒng)

本節(jié)采用基于IEEE 13 節(jié)點的交直流混合系統(tǒng)進(jìn)行測試,分別針對有載調(diào)壓器電壓比是否已知、模型的目標(biāo)函數(shù)是否為線路有功功率損耗以及平衡節(jié)點是否存在有載調(diào)壓器來改變模型,并將凸松弛模型的計算結(jié)果代入一般潮流算法中計算,對比一般潮流算法與凸松弛算法的求解結(jié)果。具體的模型調(diào)整方式如下:

(1)有載調(diào)壓器電壓比已知且凸松弛模型的目標(biāo)函數(shù)為0,將凸松弛模型中的分布式電源無功功率的計算結(jié)果代入一般潮流算法。

(2)有載調(diào)壓器電壓比未知且凸松弛模型的目標(biāo)函數(shù)為線路有功功率損耗,將凸松弛模型中的有載調(diào)壓器電壓比和分布式電源的無功功率計算結(jié)果代入一般潮流算法。

(3)有載調(diào)壓器電壓比已知且凸松弛模型的目標(biāo)函數(shù)為線路有功功率損耗,將凸松弛模型中的分布式電源無功功率的計算結(jié)果代入一般潮流算法。

(4)有載調(diào)壓器電壓比未知且凸松弛模型的目標(biāo)函數(shù)為0,將凸松弛模型中的有載調(diào)壓器電壓比和分布式電源的無功功率計算結(jié)果代入一般潮流算法。

(5)去掉平衡節(jié)點處的有載調(diào)壓器且凸松弛模型的目標(biāo)函數(shù)為0,將凸松弛模型中的分布式電源無功功率的計算結(jié)果代入一般潮流算法。

(6)去掉平衡節(jié)點處的有載調(diào)壓器且凸松弛模型的目標(biāo)函數(shù)為線路有功功率損耗,將凸松弛模型中的分布式電源無功功率的計算結(jié)果代入一般潮流算法。

將凸松弛模型的優(yōu)化結(jié)果作為已知量代入原有的一般潮流計算模型中,計算結(jié)果與原有一般潮流模型的電壓計算結(jié)果進(jìn)行比較,將原有一般潮流模型的電壓計算結(jié)果作為真實值,對比結(jié)果見表1。

表1 SOCP 凸松弛優(yōu)化計算精度對比分析Tab.1 Comparative analysis of the accuracy of SOCP relaxation optimization

由表1 可知,SOCP 凸松弛模型所得的結(jié)果與原始一般潮流模型中所解得的結(jié)果較為接近,說明求解器所計算的結(jié)果具有一定的準(zhǔn)確性與合理性,同時也證明本文中所提出的 SOCP 凸松弛模型的有效性。

為了驗證凸松弛之后所得最優(yōu)解在原約束下的等式精確性,定義原等式約束的二階錐松弛偏差矢量無窮范數(shù)為

式中,D為二階錐松弛偏差矢量無窮范數(shù);EL、ER分別為表示凸松弛前等式約束的等號左邊項與等式約束的等號右邊項。

在上述六種工況的最優(yōu)解處,有

由式(30)可以看出,本文所提出的二階錐規(guī)劃凸松弛模型是精確的,模型具有一定的合理性。

3.1.2 計及分布式電源反送功率

本節(jié)在3.1.1 節(jié)的基礎(chǔ)上進(jìn)一步分析當(dāng)分布式電源容量大于負(fù)荷時,分布式電源向電網(wǎng)反送功率對主動配電網(wǎng)潮流計算結(jié)果的影響。本節(jié)中,增加分布式電源的有功功率上、下邊界約束,下邊界設(shè)定為正常容量(不反送功率時的容量)的5 倍,上邊界設(shè)定為正常容量的6 倍,使其向電網(wǎng)反送功率,利用3.1.1 節(jié)中六種工況進(jìn)行測試,結(jié)果見表2。

表2 計及分布式電源反送功率的SOCP 凸松弛優(yōu)化計算精度對比分析Tab.2 Comparative analysis of the accuracy of SOCP relaxation optimization considering the reverse power ofdistributed generator

由表2 可知,與分布式電源不反送功率時的運(yùn)行結(jié)果相比,當(dāng)分布式電源的容量超過負(fù)荷時,分布式電源向電網(wǎng)反送功率對凸松弛優(yōu)化計算結(jié)果影響較小,說明在分布式電源出現(xiàn)反送功率的情況下,本文所提模型依然有效,進(jìn)一步證明本文所提模型的合理性與有效性。

3.1.3 對比不同求解器的求解效率

本文基于YALMIP 平臺上的商用和非商用求解器進(jìn)行簡單的求解效率對比。對比測試的電力系統(tǒng)分別為基于IEEE 13 節(jié)點的交直流混合系統(tǒng)及IEEE 123 節(jié)點系統(tǒng),本文使用的Matlab 版本為R2018b,內(nèi)存為16.0GB,CPU 為3.40GHz。本文主要比較的是各求解器對相同SOCP 凸松弛模型的運(yùn)算時間,并簡短討論求解器的求解精確度。

表3 是使用不同的求解器求解IEEE 13 節(jié)點的交直流混合系統(tǒng)凸松弛模型的求解結(jié)果,其中求解器SCS 的表現(xiàn)最好,求解時間最短,且求得的電壓幅值與實際潮流的結(jié)果相近。ECOS 的表現(xiàn)次之,求解的時間稍長,求解結(jié)果與SCS 的求解結(jié)果相近,但精度稍差。SeDuMi 和SDPT3 兩者的表現(xiàn)相近,計算時長相近,且計算結(jié)果也相似,計算精度與ECOS 相似。因為測試的算例規(guī)模較小,只有13 個節(jié)點,因此并不能看出哪個求解器有特別突出的求解能力。MOSEK 求解器無法求解。SCS 求解器高效且準(zhǔn)確的求解能力是因為它使用ADMM 分布式方法,提高了運(yùn)算效率,因此在求解大規(guī)模的二階錐規(guī)劃問題上SCS 求解器有更明顯的優(yōu)勢。

表3 不同求解器求解基于IEEE 13 節(jié)點的交直流混合系統(tǒng)時間對比Tab.3 Time comparison of different solvers for solving AC-DC hybrid system based on IEEE 13 nodes

表4 是IEEE 123 節(jié)點系統(tǒng)測試結(jié)果,比較了9種不同求解的運(yùn)算時間,其中商用求解器MOSEK的表現(xiàn)最好,其次是SDPA 求解器,但是其求解結(jié)果與實際潮流結(jié)果存在較大差異。SDPT3 和SCS 算法表現(xiàn)相似,計算時間較短,且準(zhǔn)確率較高,其他求解器表現(xiàn)一般,均存在計算誤差較大或計算時間較長等問題。

表4 不同求解器求解IEEE 123 節(jié)點系統(tǒng)時間對比Tab.4 Time comparison of different solvers for IEEE 123 nodes system

由表3 和表4 可知,針對本文IEEE 13 節(jié)點系統(tǒng)和IEEE 123 節(jié)點系統(tǒng),不同求解器的計算效率存在差異。在求解相關(guān)問題時,多數(shù)文獻(xiàn)只是固定地選擇MOSEK 求解器,并沒有具體討論不同求解器的求解效率,根據(jù)本文分析結(jié)果,說明對于不同的問題MOSEK 不一定總是表現(xiàn)最好,應(yīng)根據(jù)實際問題合理地選擇求解器。

4 結(jié)論

本文提出了考慮分布式電源、不同接線方式的有載調(diào)壓器、交直流換流器等三相模型的二階錐規(guī)劃凸松弛模型,考慮相間耦合,同時給出了分布式電源、有載調(diào)壓器等元件的物理特性,并采用多種二階錐規(guī)劃求解器求解,松弛后的優(yōu)化模型能夠保證得到全局最優(yōu)解。在求解電力系統(tǒng)最優(yōu)潮流時,可作為控制的最終控制結(jié)果,也可以作為非凸優(yōu)化求解器初值以更快得到全局最優(yōu)解,具有一定的實用性。通過IEEE 13 節(jié)點系統(tǒng)和IEEE 123 節(jié)點系統(tǒng)的計算驗證了本文模型的有效性和可行性。

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 精品综合久久久久久97超人| 亚洲成a∧人片在线观看无码| 亚洲aaa视频| 一本色道久久88亚洲综合| 国产精品久久自在自2021| 成人精品视频一区二区在线| 久久a级片| 精品国产免费观看一区| 九九热精品视频在线| 2020久久国产综合精品swag| 在线观看91香蕉国产免费| 99热这里只有精品免费| 成人日韩视频| 亚洲AV色香蕉一区二区| 永久免费无码日韩视频| 亚洲国产精品无码AV| 人妻中文久热无码丝袜| A级全黄试看30分钟小视频| 日本91视频| 无码福利日韩神码福利片| 午夜无码一区二区三区在线app| 国产午夜不卡| 国产乱子伦精品视频| 久久亚洲国产一区二区| 精品久久高清| 91小视频在线观看免费版高清| 大陆精大陆国产国语精品1024 | 97色婷婷成人综合在线观看| 亚洲男人的天堂网| av在线人妻熟妇| 亚洲无码高清视频在线观看| 亚洲成人精品| 久久国产精品影院| 国内精品伊人久久久久7777人| 国产91av在线| 亚洲综合欧美在线一区在线播放| 国产精品无码一二三视频| 亚洲av日韩综合一区尤物| 日韩无码视频播放| 日韩精品毛片人妻AV不卡| 毛片大全免费观看| 亚洲午夜国产片在线观看| 91网址在线播放| 日韩精品一区二区三区免费在线观看| 欧美高清日韩| 日本福利视频网站| 国产亚洲高清在线精品99| 精品偷拍一区二区| 亚洲精品日产精品乱码不卡| 91精品久久久无码中文字幕vr| 试看120秒男女啪啪免费| 精品久久久久久成人AV| 午夜福利在线观看入口| 国产精品久久久久鬼色| 在线观看无码av免费不卡网站 | 九色视频线上播放| 亚洲精品va| 欧美午夜在线播放| 亚洲精品自拍区在线观看| 日本在线欧美在线| 亚洲无码A视频在线| 一级成人a毛片免费播放| 国产高清国内精品福利| 伊人久久精品无码麻豆精品| 伊人精品视频免费在线| 精品国产美女福到在线不卡f| 影音先锋亚洲无码| 色综合热无码热国产| 中文字幕精品一区二区三区视频| 精品国产三级在线观看| 国产精品成人一区二区| 性色在线视频精品| 黄色污网站在线观看| 日韩AV无码免费一二三区 | 亚洲国产成人自拍| 久久精品嫩草研究院| 亚欧成人无码AV在线播放| 第一页亚洲| 国产精品网址你懂的| 欧美精品导航| 亚洲精品卡2卡3卡4卡5卡区| 国产一区自拍视频|