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

基于改進非參數核密度估計和拉丁超立方抽樣的電動公共客車負荷模型

2016-04-07 00:35:42繆鵬彬余娟史樂峰劉國平梁明李文沅任洲洋輸配電裝備及系統安全與新技術國家重點實驗室重慶大學重慶400044國網重慶市電力公司電力科學研究院重慶403
電工技術學報 2016年4期
關鍵詞:方法

繆鵬彬余 娟史樂峰劉國平梁 明李文沅任洲洋(. 輸配電裝備及系統安全與新技術國家重點實驗室(重慶大學) 重慶 400044 . 國網重慶市電力公司電力科學研究院 重慶 403)

?

基于改進非參數核密度估計和拉丁超立方抽樣的電動公共客車負荷模型

繆鵬彬1余 娟1史樂峰2劉國平2梁 明1李文沅1任洲洋1
(1. 輸配電裝備及系統安全與新技術國家重點實驗室(重慶大學) 重慶 400044 2. 國網重慶市電力公司電力科學研究院 重慶 401123)

摘要為了建立更精確的電動公共客車充電負荷模型,對起始荷電狀態和起始充電時間的概率分布和抽樣方法進行研究。首先針對概率分布,提出采用邊界核的自適應非參數核密度估計算法,該方法無需概率分布的任何假設,能有效解決邊界偏差和缺乏局部適應性的問題,從而提高了起始荷電狀態和起始充電時間概率分布的準確度和適應性。其次針對抽樣環節,提出結合三次樣條插值法的改進拉丁超立方抽樣算法,該方法不僅可以解決累積分布函數的反函數較難求解的問題,而且具有準確度高、抽樣規模小等優點。最后,以這兩種算法為基礎,建立電動公共客車充電負荷模型,并將其與實測數據和傳統參數估計方法所得結果進行對比。仿真結果證明了所提方法的適應性和有效性。

關鍵詞:電動公共客車 充電負荷 自適應非參數核密度估計算法 邊界核 拉丁超立方抽樣

國家自然科學基金(51477017),重慶市基礎與前沿研究計劃一般項目(cstc2014jcyjA90017)和國家電網公司科技項目(SGCQDK00DJJS1500056)資助。

Electric Public Bus Load Model Based on Improved Kernel Density Estimation and Latin Hypercube Sampling

Miao Pengbin1Yu Juan1Shi Lefeng2Liu Guoping2Liang Ming1Li Wenyuan1Ren Zhouyang1(1. State Key Laboratory of Power Transmission Equipment & System Security and New Technology
Chongqing University Chongqing 400044 China
2. Chongqing Electric Power Research Institute Chongqing 401123 China)

Abstract To build a more precise charging load model of electric public buses, the probability distribution and sampling method of charging state and charging start time are studied. Firstly, an adaptive kernel density estimation with boundary Kernel algorithm is proposed to build probability distribution models. The proposed algorithm does not require any assumptions about the probability distribution, and can also solve the problems of boundary bias and lacking of local adaptability, which improves the precision and adaptability of probability distribution. Secondly, Latin hypercube sampling with cubic spline interpolation algorithm is proposed to figure out the inverse cumulative distribution function of probability distribution. The proposed algorithm has high precision and sampling efficiency. Finally, based on these two algorithms, electric public bus charging load model is set up. The simulation result demonstrates the effectiveness and adaptability of the proposed method.

Keywords:Electric public bus, charging load, adaptive Kernel density estimation, boundary Kernel, Latin hypercube sampling

0 引言

電動汽車作為新一代的交通工具,具有低能耗、低污染的巨大優勢,在可以預見的未來將迎來飛速的發展。對電動汽車的充電負荷進行建模是研究電動汽車對電網影響、充電基礎設施規劃以及電網規劃與運行的前提[1,2]。電動汽車目前主要分為私家車、公共客車和出租車三類,其中電動公共客車的發展較為快速。電動公共客車一般集中在大型充電站進行快速充電,配套有較為完善的充電數據采集系統,能夠為基于概率統計方法建立其充電負荷模型提供充足的歷史數據。

目前電動汽車充電負荷模型研究廣泛采用基于蒙特卡洛模擬的統計學建模方法[3-6]。該方法是在建立諸如起始荷電狀態(State Of Charge, SOC)、起始充電時間等概率分布的基礎上,考慮充電功率、電動汽車數量及電池容量等影響因素,運用蒙特卡洛數學試驗累加每輛電動汽車的充電負荷得到總體充電負荷。考慮到電動公共客車一般集中在大型充電站進行快速充電(充電功率恒定,充至滿電),且具有電池容量標配,日充電次數基本恒定等特點,因此起始SOC和起始充電時間的精確概率分布是決定電動公共客車充電負荷模型的關鍵因素。

目前一般采用傳統參數估計方法建立起始SOC和起始充電時間的概率分布,即預先假設其概率密度函數服從某種分布,然后依據歷史樣本數據來估計該分布的參數。對于起始SOC,大部分文獻假設其服從正態分布[3,7],文獻[5]則假設其服從廣義極值分布;對于起始充電時間,大部分文獻假設其服從均勻分布[6,8],文獻[3]則根據不同車輛類型,假定其服從均勻分布、正態分布和直角梯形分布,文獻[5]假定其服從威布爾分布。這些概率分布模型雖有不同,但是均屬于傳統參數估計范疇。如果假設的概率分布與實際情況不符,那么必然會導致計算所得的充電負荷模型偏差較大。因此,引入非參數核密度估計(Kernel Density Estimation, KDE)算法是十分必要的,該算法不需要概率分布的任何假設,能夠揭示隱藏在歷史樣本數據中的統計信息,有效提高概率分布準確度。需要指出的是,本文所分析的起始SOC和起始充電時間都是有上、下界約束的數據,利用KDE算法求解其概率分布時會在上、下邊界處出現邊界偏差(隨機變量上、下邊界以外的概率密度值仍然大于0,從而造成誤差),這可以通過采用邊界核的KDE算法來解決[9];此外,數據分布的光滑程度不具有一致性,KDE算法假設樣本區間內光滑程度一致(采用固定帶寬)較為理想化,與實際情形不符,對此B. W. Silverman提出了自適應核密度估計法[10],從而解決了KDE算法存在的缺乏局部適應性問題。相較傳統參數估計方法,采用自適應核密度估計法對于不同充電場景下的充電負荷特性具有更好的適應性。本文將邊界核與自適應核密度估計法相結合,提出基于邊界核的自適應非參數核密度估計(Adaptive Kernel Density Estimation with Boundary,AKDEB)算法,以有效解決KDE算法存在的邊界偏差和缺乏局部適應性的問題。

基于非參數核密度估計算法的概率密度函數是求和形式,其累積分布的反函數難以求解。因此,在非參數核密度估計抽樣中,不依賴于累積分布反函數的舍選法得到了一定的應用[11],而相比舍選法抽樣規模更小、準確度更高的拉丁超立方抽樣(Latin Hypercube Sampling,LHS)方法則無法直接應用[12]。然而,基于一組已知的樣本值和對應的累積概率,利用插值法便可以求解出任意累積概率所對應的樣本值,三次樣條插值法由于計算簡便、光滑性好和準確度高等優點得到了廣泛應用[13]。因此,本文提出一種基于三次樣條插值法的改進拉丁超立方抽樣(LHS with Cubic Spline Interpolation,LHS-CSI)法,以解決累積分布反函數難以求解時不能應用LHS方法的問題,并且保留了LHS方法抽樣規模小、準確度高的優點。

綜上,本文首先提出AKDEB算法建立起始SOC和起始充電時間的概率分布模型;其次,提出LHS-CSI算法對起始SOC和起始充電時間進行抽樣,并建立電動公共客車充電負荷模型;最后,以重慶市渝北空港電動客車專用充電站的實際統計數據為基礎進行仿真分析,以驗證本文所提方法的適應性和有效性。

1 基于AKDEB估計算法的起始SOC和起始充電時間概率分布

本文所提的AKDEB算法結合了邊界核和自適應核密度估計的優點,可以很好地解決KDE算法存在的邊界偏差和缺乏局部適應性的問題。

下面介紹該算法的具體原理:假設X1, X2,…, Xn為n個樣本,定義域為[a,b],概率密度函數為f(x),則此概率密度函數的KDE估計[8]為

常規KDE算法存在邊界偏差問題,可以通過使用“邊界核”這一特殊的核函數來修正[9]。邊界核的表達式為

式中

將式(1)中的K(Z)用式(2)表示的邊界核B(Z)進行替換,即得到采用邊界核的核密度估計為

采用式(3)的核密度估計僅僅能解決邊界偏差,但是由于采用的是固定帶寬h,因此并不能根據局部區域的樣本數據調節帶寬,從而存在缺乏局部適應性的問題,對此B. W. Silverman提出了自適應核密度估計法[10],即把式(1)中的h用式(4)代替。式中,f

(Xi)一般采用式(1)的估計值作為先導估計。hv則是由交叉驗證法[10]得到的最優帶寬。

本文將邊界核與自適應核密度估計算法相結合,提出AKDEB算法,即把式(3)中的固定帶寬h用式(4)的可變帶寬h(Xi)代替,則概率密度函數的AKDEB估計為

需要注意的是,由于式(5)中的邊界核結合了自適應核密度估計算法, 因此原來的相應地轉變為。

以起始SOC和起始充電時間的歷史樣本數據為基礎,利用式(5)便可以得到其概率密度估計。起始SOC和起始充電時間在邊界以外的概率密度估計值為0,在邊界以內則會根據局部區域的樣本數據調節帶寬,這樣就可以解決常規KDE算法的邊界偏差和缺乏局部適應性問題。

2 基于LHS-CSI抽樣算法的電動公共客車充電負荷模型

首先對本文提出的LHS-CSI算法的基本原理進行介紹;然后,基于LHS-CSI抽樣算法對起始SOC和起始充電時間進行隨機抽樣,建立電動公共客車充電負荷模型。

2.1 LHS-CSI抽樣算法

對于某一隨機變量X,定義域為[a,b],累積概率分布函數為

設N代表采樣規模,采用中值拉丁超立方抽樣法將曲線Y=F(X)的縱軸分成N個等間距不重疊的區間,選擇每個區間的中點作為Y的采樣值yk=(k-0.5)/N (k=1,2,…,N),然后用函數Y=F(X)的反函數來計算隨機變量X的第k個采樣值為

以上介紹的是常規LHS的采樣方法,但是本文中F-1(yk)的顯式表達式難以求解,因此需要借助三次樣條插值法[14]:

(1)將隨機變量X的定義域[a,b]均分為m個點(a=x1<x2<…<xm-1<xm=b),并根據式(5)求出點xi(i=1,2,…,m)處的概率密度估計值,從而點xi處的概率估計值近似等于因此點xi處的累加概率估計值為

(3)對于任一累積概率值yk=(k-0.5)/N,k= 1,2,…,N,一定落在某個小區間內,將yk代入式(9)則可以求解出樣本值。

經過以上步驟,便可以得到N個起始SOC和起始充電時間的隨機樣本,為建立電動公共客車充電負荷模型提供數據輸入。LHS-CSI算法準確度高、抽樣效率高的特點可以保證最后建立電動公共客車充電負荷模型的準確性和高效性。

2.2 電動公共客車充電負荷模型

假設充電站一天內接納M車次(同一輛電動公共客車一天內可能充電多次,M車次表示充電站一天內所有電動公共客車的充電總次數)的待充車輛,所有待充電動公共客車的電池容量為Cap,充電過程近似為恒功率充電,充電功率為Pc。對于任意第i (i=1,2,…,M)車次的待充公共客車,由2.1節可以得到它對應的起始SOC和起始充電時間樣本值,記作(Si,ti)。根據Si可以直接計算出其充電持續時間為

根據式(10)得到的充電持續時間Ti,可以得出該車輛的充電時間區間為[ti, ti+Ti]。定義電動公共客車的上班時間為tstart,下班時間為tend,則上班時間內任一時刻t的充電負荷可以表示為

式中

利用式(11)可以求解出任意時刻t的充電負荷,當求解出區間[tstart, tend]內所有時刻點的充電負荷并連接成線時,即可求出最終的充電日負荷曲線。實際應用中,一般以Δt(Δt可以取為1s、1min和1h等)為時間間隔將區間[tstart, tend]分為若干個小區間,求解出小區間上的充電負荷并連接成線,形成最終的日負荷曲線,具體計算流程如圖1所示。

3 算例分析

3.1 數據描述

本文采用重慶市渝北空港電動客車專用充電站的實際統計數據(2013年3月~2014年4月)來證明所提方法。該充電站電動公共客車采用Pc=260kW恒功率充電。電池容量Cap=100kW·h。上班時間tstart=420min(上午7點),下班時間tend=1 290min(晚上9點半)。起始SOC的定義域為[0,100],m=101(起始SOC概率密度直方圖的區間間隔為1%)。起始充電時間的定義域為[420min,1 290min],m=871(起始充電時間概率密度直方圖的區間間隔為1min)。充電站日接納車次M=100。采樣規模N=10 000。充電負荷曲線計算時間步長為Δt =1min。Matlab運行的計算機環境為Inter Core 2 Duo核處理器,內存8G。

圖1 電動公共客車充電負荷計算流程Fig.1 Flowchart of the electric public bus charging load

3.2 起始SOC和起始充電時間概率分布的仿真結果分析

本文采用卡方檢驗和后驗檢驗[15]作為評估概率分布擬合優度的指標。表1和表2給出了起始SOC和起始充電時間的常用概率分布的卡方檢驗結果。卡方統計量越小表明概率分布的擬合效果越好。從表中可以看出,在幾種常用的傳統參數估計方法中,正態分布對于起始SOC的擬合相對較好,均勻分布對于起始充電時間的擬合誤差也在可接受的范圍內,因此這兩種分布模型的廣泛使用具有一定的合理性。然而,采用傳統參數分布時其估計結果均不能通過嚴格統計意義的卡方檢驗;相比之下,采用非參數核密度估計方法的卡方統計量均遠小于傳統參數估計方法的統計量,而且能通過卡方檢驗。本文提出的AKDEB算法比KDE算法的卡方統計量更小,其擬合效果更優。

表1 起始SOC四種概率分布的卡方檢驗結果Tab.1 Chi square test results of four probability distribution models with SOC

表2 起始充電時間三種概率分布的卡方檢驗結果Tab.2 Chi square test results of three probabilitydistribution models with charging start time

圖2和圖3分別給出了起始SOC和起始充電時間的概率密度曲線和直方圖。從圖2和圖3的后驗檢驗也可以看出本文所提的概率分布較之傳統參數分布具有更好的擬合準確度。綜上,本文提出的AKDEB算法可以有效提高概率分布模型的準確度。

由圖2和圖3還可以得出以下兩個結論:

(1)常規KDE算法形成的概率密度曲線確實存在邊界偏差,起始SOC<0和SOC>100%處的概率密度值仍然是>0的,起始充電時間在<420min和>1 290min處的概率密度值也是>0的,而這些都是與實際情況不相符的,相比之下,本文所提出的AKDEB算法修正了邊界偏差,符合實際情況,進一步提高了擬合準確度。

圖2 起始SOC概率密度曲線和直方圖Fig.2 The probability density curve and histogram of SOC

圖3 起始充電時間概率密度曲線和直方圖Fig.3 The probability density curve and histogram of charging start time

(2)從圖3可以看出AKDEB算法所得的概率密度曲線中部較為粗糙,即曲線中部的帶寬值較小,而曲線兩端則較為光滑,即曲線兩端的帶寬值較大,從而實現了帶寬值根據樣本數據的變化進行調節,解決了KDE算法缺乏局部適應性的問題。

3.3 起始SOC和起始充電時間抽樣的仿真結果分析

圖4 LHS-CSI算法和舍選法準確度對比Fig.4 Accuracy comparison between LHS-CSI and acceptance-rejection sampling algorithm

本文以起始SOC為例,采用方差系數作為衡量LHS-CSI算法和舍選法準確度的指標。圖4、圖5給出了本文抽樣方法和舍選法在準確度和抽樣時間上的對比。在抽樣準確度方面,如圖4所示,如果以方差系數1%作為收斂判據,那么LHS-CSI算法在抽樣次數為500次左右時即可收斂,而舍選法在抽樣次數達到5 000次時仍未收斂,經過進一步計算,抽樣次數達到20 000次時舍選法才能收斂。如果以方差系數5%作為收斂判據,那么LHS-CSI算法在抽樣次數為20次時即可收斂,而舍選法在抽樣次數達到1 500次時才能收斂。在抽樣時間方面,如圖5所示,由于LHS-CSI算法的抽樣時間主要集中在計算典型點xi(a=x1<x2<…<xm-1<xm=b)處的概率密度估計值,因此在m選定的情況下抽樣時間基本保持不變。但是,舍選法每次抽樣都需要進行概率密度值的計算,因此其計算時間會隨著抽樣次數的增加呈線性增長的趨勢。綜上,LHS-CSI算法在準確度和抽樣時間上均優于舍選法。

圖5 LHS-CSI算法和舍選法抽樣時間對比Fig.5 Time comparison between LHS-CSI and acceptance-rejection sampling algorithm

3.4 電動公共客車充電日負荷曲線的仿真結果分析

圖6 基于AKDEB算法與傳統參數估計算法的電動公共客車充電日負荷曲線Fig.6 Daily load curves of the electric public buses based on AKDEB and traditional parameter estimation algorithm

圖6給出了三條充電日負荷曲線,依次對應于實測數據、本文所提方法、傳統參數估計方法,其中,傳統參數估計方法中對起始SOC和起始充電時間分別采用的是正態分布和均勻分布。從三條日負荷曲線的對比可以看出,采用本文方法對電動汽車充電負荷進行建模時與實際情況較為接近,其曲線走勢基本與實測數據相同,而傳統參數統計方法仿真結果與實測數據之間的誤差較大,其誤差最大的兩個區域集中在上午9∶20和晚上8∶40左右。

表3給出了以平均相對誤差、峰(谷)時相對誤差為指標的兩種方法的準確度對比。從表中可以看出,采用傳統參數估計方法所得的充電日負荷曲線與實測數據的誤差均遠大于本文所提算法的誤差,傳統參數估計方法的谷時相對誤差更是達到了208.9%。綜上,本文所提方法可以有效提高電動公共客車充電負荷模型的準確度。

表3 本文方法和傳統參數估計方法的準確度對比Tab.3 Accuracy comparison between AKDEB and traditional parameter estimation algorithm

4 結論

首先,本文提出了AKDEB估計算法,并利用該算法對起始SOC和起始充電時間的概率分布進行了研究。該算法不需要概率分布的任何假設,能夠揭示隱藏在歷史樣本數據中的統計信息,有效提高了概率分布的準確度。此外,該算法將邊界核和自適應核密度估計方法的優點有效結合,解決了常規KDE算法存在的邊界偏差和缺乏局部適應性的問題。其次,本文提出了LHS-CSI抽樣算法,對起始SOC和起始充電時間進行隨機抽樣。該算法解決了累積分布函數的反函數難以求解顯示表達式時不能直接應用LHS抽樣方法的問題,且有效保留了LHS方法抽樣準確度高、計算時間短的優點。最后,基于本文所提的以上兩種算法,建立了電動公共客車充電負荷模型。本文所提出方法準確度高、適應性強且效率高,保證了最終所得到的充電負荷模型的準確性、適應性和高效性,為研究電動汽車對電網影響、充電基礎設施規劃以及電網規劃與運行等問題打下了良好的基礎。仿真結果證明了本文所提方法的適應性和有效性。

參考文獻

[1] 葛文捷, 黃梅, 張維戈. 電動汽車充電站經濟運行分析[J]. 電工技術學報, 2013, 28(2): 15-21.

Ge Wenjie, Huang Mei, Zhang Weige. Economicoperation analysis of the electric vehicle charging station[J]. Transactions of China Electrotechnical Society, 2013, 28(2): 15-21.

[2] 肖湘寧, 溫劍鋒, 陶順,等. 電動汽車充電基礎設施規劃中若干關鍵問題的研究與建議[J]. 電工技術學報, 2014, 29(8): 1-10.

Xiao Xiangning, Wen Jianfeng, Tao Shun, et al. Study and recommendations of the key issues in planning of electric vehicles’ charging facilities[J]. Transactions of China Electrotechnical Society, 2014, 29(8): 1-10.

[3] 羅卓偉, 胡澤春, 宋永華, 等. 電動汽車充電負荷計算方法[J]. 電力系統自動化, 2011, 35(14): 36-42.

Luo Zhuowei, Hu Zechun, Song Yonghua,et al. Study on plug-in electric vehicles charging load calculating[J]. Automation of Electric Power Systems, 2011, 35(14): 36-42.

[4] 楊冰, 王麗芳, 廖承林. 大規模電動汽車充電需求及影響因素[J]. 電工技術學報, 2013, 28(2): 22-27.

Yang Bing, Wang Lifang, Liao Chenglin. Research on power-charging demand of large-scale electric vehicles and its impacting factors[J]. Transactions of China Electrotechnical Society, 2013, 28(2): 22-27. [5] Pashajavid E, Golkar M A. Charging of plug-in electric vehicles: Stochastic modelling of load demand within domestic grids[C]//2012 20th Iranian Conference on Electrical Engineering, 2012: 535-539.

[6] 徐國鈞, 劉永勝, 李題印, 等. 基于層次分析和概率模擬的電動汽車對配網負荷影響研究[J]. 電力系統保護與控制, 2012, 40(22): 38-45.

Xu Guojun, Liu Yongsheng, Li Tiyin, et al. Study on the impact of electric vehicles charging load integrated into distribution network based on analytic hierarchy process and probabilistic simulation[J]. Power System Protection and Control, 2012, 40(22): 38-45.

[7] Leou Rong Ceng, Su Chun Lien, Lu Chan Nan. Stochastic analyses of electric vehicle charging impacts on distribution network[J]. IEEE Transactions on Power Systems, 2014, 29(3): 1055-1063.

[8] Dai Qian, Cai Tao, Duan Shanxu, et al. Stochastic modeling and forecasting of load demand for electric bus battery-swap station[J]. IEEE Transactions on Power Delivery, 2014, 29(4): 1909-1917.

[9] Jeffery S Simonoff. Smoothing methods in statistics[M]. New York: Springer-Verlag, 1996.

[10] Silverman B W. Density estimation for statistics and data analysis[M]. Boca Raton: CRC Press, 1986.

[11] 趙淵, 張夏菲, 周家啟. 電網可靠性評估的非參數多變量核密度估計負荷模型研究[J]. 中國電機工程學報, 2009, 29(31): 27-33.

Zhao Yuan, Zhang Xiafei, Zhou Jiaqi. Load modeling utilizing nonparametric and multivariate Kernel density estimation in bulk power system reliability evaluation[J]. Proceedings of the CSEE, 2009, 29(31): 27-33.

[12] Mckay M D, Beckman R J, Conover W J. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J]. Technimetrics, 1979, 21(2): 239-245.

[13] Young T, Mohlenkamp M J. Introduction to numerical methods and Matlab programming for engineers[M]. Athens: Ohio University, 2012.

[14] 張德豐, 等. MATALAB 數值計算方法[M]. 北京:機械工業出版社, 2010.

[15] Qin Zhilong, Li Wenyuan, Xiong Xiaofu. Estimating wind speed probability distribution using Kernel density method[J]. Electric Power Systems Research, 2011, 81(12): 2139-2146.

繆鵬彬 男,1991年生,碩士研究生,研究方向為電動汽車與電網安全穩定分析。

E-mail: miaopengbin@qq.com

余 娟 女,1980年生,博士,副教授,研究方向為無功優化、風險評估、電壓穩定與電力系統靜態等值。

E-mail: cqu.juan.yu@gmail.com(通信作者)

作者簡介

收稿日期2014-09-23 改稿日期 2015-09-17

中圖分類號:TM743; U469.72+1

猜你喜歡
方法
中醫特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數學教學改革的方法
河北畫報(2021年2期)2021-05-25 02:07:46
化學反應多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 干中文字幕| 亚洲性影院| 日韩中文字幕免费在线观看| 又粗又大又爽又紧免费视频| 国产三级毛片| 日韩精品一区二区三区中文无码| 久草性视频| 亚洲综合色婷婷| 91无码人妻精品一区| 在线五月婷婷| www精品久久| 欧美日韩国产系列在线观看| 亚洲国产日韩一区| 色综合成人| 免费在线观看av| 国产福利大秀91| 久久免费精品琪琪| 无码丝袜人妻| 人妻一区二区三区无码精品一区| 潮喷在线无码白浆| 毛片网站在线看| 又爽又大又光又色的午夜视频| 国产乱子伦精品视频| 亚洲一区二区黄色| 九九九精品视频| 伊人久久影视| 青青青国产视频手机| 日韩国产高清无码| 国产91九色在线播放| 国产日韩欧美精品区性色| 67194亚洲无码| 久久一色本道亚洲| 在线va视频| 欧美激情视频二区| 国产又黄又硬又粗| 一区二区三区国产| 全部毛片免费看| 欧美亚洲国产精品第一页| 午夜精品一区二区蜜桃| 国产爽妇精品| 国产微拍一区| 国产日产欧美精品| 人人爽人人爽人人片| 日本少妇又色又爽又高潮| 尤物亚洲最大AV无码网站| 国产黄色视频综合| 久久国产精品国产自线拍| 国产日韩av在线播放| 老汉色老汉首页a亚洲| 伊人成人在线| 欧美亚洲欧美区| 精品伊人久久大香线蕉网站| 亚洲国产精品无码久久一线| 精品91视频| 亚洲第一黄色网址| 亚洲精品成人片在线播放| 欧美精品在线视频观看| 亚洲国产成人久久精品软件 | 香蕉视频在线观看www| 久久精品免费国产大片| 亚洲熟女中文字幕男人总站| 无码精品一区二区久久久| 青青草原国产精品啪啪视频| 国产乱人伦偷精品视频AAA| 色综合天天综合中文网| 久操中文在线| 日韩天堂网| 久久久国产精品免费视频| 手机成人午夜在线视频| 色婷婷电影网| 久久精品电影| 成人午夜天| 全部毛片免费看| 免费精品一区二区h| 久久成人免费| 国产无人区一区二区三区| 天堂网亚洲系列亚洲系列| 国产女人综合久久精品视| 亚洲婷婷在线视频| 欧洲欧美人成免费全部视频| 天天干伊人| 国产成人高清精品免费|