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

利用最小二乘與最優(yōu)控制的水輪機模型參數(shù)辨識

2020-12-01 03:15:00侯睿曾云王偉李敏古志吳一凡
軟件導刊 2020年10期

侯睿 曾云 王偉 李敏 古志 吳一凡

摘 要:為了改善一般常用水輪機模型基于實測數(shù)據(jù)參數(shù)進行辨識過程中底層編程繁瑣的缺點,基于電力系統(tǒng)仿真中常用的理想水輪機數(shù)學模型,取134組模擬實測輸入輸出數(shù)據(jù),用最小二乘與基于數(shù)值最優(yōu)化技術(shù)的最優(yōu)控制方法分別進行辨識從而獲取參數(shù),并利用兩種方法對同一實測數(shù)據(jù)辨識結(jié)果進行仿真對比。算例仿真結(jié)果表明,在給定階躍擾動下,利用最小二乘與最優(yōu)控制方法進行參數(shù)辨識后的模型輸出響應與實測響應結(jié)果(基于模擬樣本數(shù)據(jù))誤差均低于0.04,且基于數(shù)值最優(yōu)化技術(shù)的最優(yōu)控制辨識效果與實測數(shù)據(jù)擬合程度更精確,驗證了基于數(shù)值最優(yōu)化技術(shù)的最優(yōu)控制方法在線性水輪機模型參數(shù)辨識上應用的可行性,為水輪機模型參數(shù)辨識方法提供了一種可提升辨識效率的人機對話參考。

關(guān)鍵詞:最小二乘算法;線性水輪機模型;MATLAB;最優(yōu)控制方法;參數(shù)辨識

DOI:10. 11907/rjdk. 201173

中圖分類號:TP301文獻標識碼:A 文章編號:1672-7800(2020)010-0088-06

Abstract:In order to improve the disadvantages of the low level programming in the process of parameter identification based on the measured data of common hydraulic turbine models,based on the ideal hydraulic turbine mathematical model commonly used in power system simulation, 134 groups of input and output data of simulation measurement are used to identify the parameters with the least square method and the optimal control method based on numerical optimization technology, and the two methods are used to compare the simulation results of the same measured data after identification. The simulation results show that under the given step disturbance, the error between the model output response and the measured response (based on the simulated sample data) after the parameter identification of the least square and the optimal control method is less than 0.04, and the identification effect of the optimal control based on the numerical optimization technology and the measured data fit more accurately. The feasibility of the application of the optimal control method based on the numerical optimization technology in the parameter identification of the linear hydraulic turbine model is verified, which provides a man-machine conversation reference to improve the identification efficiency for the parameter identification method of the hydraulic turbine model.

Key Words:least-square algorithm;linear hydraulic turbine model;MATLAB;method of optimal controlling;parameter identification

0 引言

水力發(fā)電在我國國民發(fā)電總量中的占比越來越重,水輪機作為水力發(fā)電主要設備之一,相關(guān)研究一直備受關(guān)注。由于難以在水電站進行現(xiàn)場實驗,進行普通真機實驗可能性較低,故通過水輪機模型仿真試驗成為主要研究方式之一,建立水輪機這一被控對象的精確化模型成為研究熱點。盡管現(xiàn)有水輪機模型研究已相當成熟,且對水輪機模型參數(shù)選擇已有詳盡歸納[1-2],但由于水輪機模型參數(shù)一般是理想化計算,所以盡管水輪機模型可滿足普通仿真實驗,卻并不能較好地貼切水電站實測響應,因此用理想化計算參數(shù)模型進行電力系統(tǒng)等仿真實驗難以避免較大誤差[3-4]。故本文提供一種不同于解析建模而根據(jù)現(xiàn)場實測數(shù)據(jù)(基于模擬樣本數(shù)據(jù))辨識水輪機模型參數(shù)的方法。將水輪機及其調(diào)節(jié)系統(tǒng)統(tǒng)一進行辨識一并得到優(yōu)化的控制參數(shù),根據(jù)實測數(shù)據(jù),對水輪機及引水系統(tǒng)分環(huán)節(jié)辨識、精細化建模成為水輪機及其調(diào)節(jié)系統(tǒng)自適應優(yōu)化控制的數(shù)學模型基石[5-6]。通常系統(tǒng)辨識有4種主要目的:①估計具有特定物理意義、可表征系統(tǒng)行為的重要參數(shù);②建立可模仿真實系統(tǒng)行為的模型,使仿真能夠反映真實的系統(tǒng)特性;③以目前可測輸入輸出預測未來演變;④辨識得到描述系統(tǒng)動態(tài)特性的數(shù)學模型以利于控制器設計[7]。其中對于水輪機,其水力動態(tài)特性主要表現(xiàn)在引水系統(tǒng)動態(tài)特性中,具有明確物理意義的水流慣性時間常數(shù)是惡化機組調(diào)節(jié)品質(zhì)的關(guān)鍵因素,但通常未對該參數(shù)進行實測,一般通過理論近似計算得到,有一定誤差。辨識得到符合真實系統(tǒng)且較小水流慣性的時間常數(shù)模型,對于改善調(diào)節(jié)品質(zhì)具有一定意義[8]。基于粒子群、引力搜索、差分進化等具有強大的非線性尋優(yōu)能力的智能算法在模型參數(shù)辨識中可發(fā)揮良好作用,因此將參數(shù)辨識問題轉(zhuǎn)化為參數(shù)優(yōu)化問題已成為趨勢[9]。

研究表明,對于線性系統(tǒng)參數(shù)辨識采用應用廣泛且非常成熟的最小二乘法亦可得到良好的辨識結(jié)果。隨著MATLAB這一強大的數(shù)值仿真計算機工具的出現(xiàn),對于復雜程度不高的參數(shù)辨識問題不必基于底層算法編程才可得以實現(xiàn),有效利用MATLAB及其Simulink工具,設計人機交互界面,也可在參數(shù)辨識等工作中獲得良好的辨識結(jié)果,并形成具有一定普適性的人機交互辨識對話框[10-15]。文獻[13]基于OCD程序?qū)Ρ豢貙ο罂刂谱髁舜罅垦芯浚珣糜谒啓C模型參數(shù)辨識的研究較少;文獻[16]將最小二乘法應用于轉(zhuǎn)向系統(tǒng)的參數(shù)辨識,取得了良好效果;文獻[17-21]成功地將最小二乘法應用于水輪機模型傳遞系數(shù)求取及曲面擬合等。

綜上所述,本文基于最小二乘法成熟且操作簡單的優(yōu)秀特性,對水輪機線性系統(tǒng)參數(shù)進行辨識,進一步參考控制領(lǐng)域?qū)<已Χㄓ罱淌赱13]總結(jié)并提出的基于跟蹤誤差指標的最優(yōu)控制器設計程序(optimal controller designer),對其Simulink模塊環(huán)節(jié)進行改進以適應于模型參數(shù)辨識工作,發(fā)現(xiàn)對于根據(jù)實測輸入輸出數(shù)據(jù)進行模型參數(shù)辨識亦可以得到良好的效果,與最小二乘法參數(shù)辨識仿真結(jié)果比較,驗證其與實測數(shù)據(jù)擬合度。

1 水輪機模型及參數(shù)辨識

1.1 水輪機模型

1992年,IEEE Working Group提出了一種廣泛應用于電力系統(tǒng)分析的簡化的線性化模型。

非理想水輪機流量和出力可表示為:

其中,[Δpt]是水輪機功率增量相對值,[Δq]是水輪機流量增量相對值,[Δy]是主接力器增量相對值,系數(shù)aij是傳遞系數(shù)。

在機組并網(wǎng)運行條件下,速度偏差很小,尤其是機組并入大電網(wǎng)運行時。因此,可以忽略Δx的影響。則有:

該理想水輪機模型具有實際物理意義的參數(shù)較少,用于參數(shù)辨識的實例研究可使問題簡化。

1.2 參數(shù)辨識

辨識的目的是將實驗得到的輸入輸出數(shù)據(jù),根據(jù)目標函數(shù),從紛復繁雜的數(shù)學模型中確定出數(shù)學模型結(jié)構(gòu)。隨著模型的參數(shù)化,只有當參數(shù)均被估計出來,模型才最終被建立,該問題為參數(shù)估計。為簡化辨識問題,在單輸入單輸出系統(tǒng)(SISO系統(tǒng))中根據(jù)輸入輸出的數(shù)據(jù)序列和已知的模型結(jié)構(gòu),結(jié)合估計理論,選擇較優(yōu)、較簡明的最小二乘算法,然后進行參數(shù)估計運算。

為了獲得對象參數(shù)[θ]的估計值[θ],一般采用逐步逼近的方法。[y(k)]為系統(tǒng)實際輸出值,[u(k)]為系統(tǒng)實際輸入值。在[k]時刻,為得到此時系統(tǒng)輸出預報值[y(k)],需根據(jù)[k]的前一時刻參數(shù)估計值[θ(k-1)]及當前時刻系統(tǒng)輸入輸出矩陣值[?(k)]。其計算公式為:[y(k)=?T(k)θ(k-1)]。此時預報誤差為[e(k)=y(k)-y(k)]。將預報誤差[e(k)]引入辨識的算法中,并且根據(jù)某準則條件得出此時參數(shù)估計值[θ(k)],將其作為下一時刻參數(shù)估計值[θ]并參與到下一時刻的系統(tǒng)輸出預報值[y(k)]計算中。依次不斷循環(huán)迭代至準則函數(shù)取最小值時為止。此時,最終系統(tǒng)輸出預報值[y(∞)]和系統(tǒng)真實輸出值[y(∞)]最為逼近,獲得此時參數(shù)估計值[θ(∞)]作為最終參數(shù)辨識結(jié)果。

其參數(shù)辨識原理如圖1所示。

2 最小二乘法與最優(yōu)控制方法

2.1 最小二乘算法辨識基礎(chǔ)

大約在1795年,高斯在其著名的星體運動軌道預報研究工作中提出了最小二乘法。由于最小二乘法原理簡明、易于理解、易于編程,且收斂較快,逐漸成為估計理論的奠基石,在參數(shù)估計中應用甚廣。但在水輪機模型參數(shù)辨識問題上,最小二乘法的應用不多。本部分簡要介紹最小二乘辨識算法基本原理。

考慮上述擴展自回歸模型(ARX),若系統(tǒng)模型是傳遞函數(shù)形式,將其轉(zhuǎn)化為形如ARX模型的形式,再將其轉(zhuǎn)換為最小二乘形式。對于單輸入單輸出(SISO)的數(shù)學模型:

一般辨識既要確定模型結(jié)構(gòu),即確定階次n,也要確定參數(shù)[ai,bi],即辨識出參數(shù)。本文討論結(jié)構(gòu)參數(shù)n和d均已知的情況,即模型結(jié)構(gòu)已知情況下僅對參數(shù)進行辨識。

將式(12)寫成最小二乘形式。

其中,β為測量總次數(shù),一般來說測量總次數(shù)β應遠大于待估參數(shù)的數(shù)量。當J為最小值時的參數(shù)估計值,[θ]即為最小二乘參數(shù)辨識結(jié)果。要求得J的最小值,即為求極值問題,使J對[θ]的求導為0。

式(17)即為所要構(gòu)造的目標函數(shù)。當[Φ]滿秩時,觀測矩陣[Φ]不是一個方陣,但[ΦTΦ]是一個方陣,[Φ]滿秩,[ΦTΦ]為可逆矩陣,故容易求得:

而且此時的J對[θ]的二階導數(shù)大于0,J對[θ]的二階導數(shù)為:

故[θ]正是最終要求辨識出的參數(shù)值。

其中,Y為β組輸出的列向量矩陣形式,[Φ]為輸入輸出矩陣φ的列向量矩陣形式。

本文采用的是離線辨識方法,即采集到所有需要的輸入輸出數(shù)據(jù)(基于模擬樣本數(shù)據(jù))后對數(shù)據(jù)用批處理最小二乘法進行集中處理。

2.2 最優(yōu)控制方法在參數(shù)辨識中的應用及流程

最優(yōu)控制器在優(yōu)化控制系統(tǒng)控制參數(shù)上有較廣泛的應用,參數(shù)辨識問題屬于控制領(lǐng)域研究范疇,故將其應用到模型參數(shù)辨識問題上是有依據(jù)的。

參數(shù)辨識問題和控制問題在實質(zhì)上均可以轉(zhuǎn)換成尋優(yōu)問題,即目標誤差函數(shù)取最小值時求得需要的參數(shù)值。基于數(shù)值最優(yōu)化技術(shù)的最優(yōu)控制器設計方法不拘泥于傳統(tǒng)最優(yōu)控制格式,目標函數(shù)可以任意定義,因此應用前景更好。參考薛定宇教授[13]提出的基于跟蹤誤差指標的最優(yōu)控制器設計程序,用MATLAB與其中的Simulink模塊求得參數(shù)最優(yōu)辨識值。用Simulink搭建待辨識模型與原真實模型的并聯(lián)系統(tǒng),附加原系統(tǒng)模型給定的階躍擾動和噪聲作為輸入模擬真實系統(tǒng),將附加待辨識模型給定的階躍擾動作為輸入。其中,待辨識模型允許將待辨識參數(shù)設為未知數(shù)a1、a2等格式。將待辨識系統(tǒng)的輸出與原真實系統(tǒng)模型的輸出連接,并求誤差給ITAE模塊。ITAE準則即為時間乘以誤差絕對值積分的誤差積分準則,需要估計ITAE指標收斂時間并將其作為參數(shù)辨識過程的仿真時間。理論上應該選擇正無窮的時間作為仿真時間,但顯然不符合實際情況。經(jīng)研究表明,選擇ITAE曲線進入穩(wěn)態(tài)后的時間1~2倍以內(nèi)均可滿足仿真要求,對待辨識參數(shù)值的影響不大,且時間選取過長將影響暫態(tài)結(jié)果。

在仿真過程中,Simulink模塊搭建完成后,需調(diào)用OCD最優(yōu)參數(shù)辨識對話界面程序,將Simulink模塊名輸入模型欄,將待辨識參數(shù)名輸入至參數(shù)欄,可根據(jù)人為約束設定待辨識參數(shù)上下界,設定仿真時間,在仿真時間內(nèi),OCD程序會不斷改變待辨識參數(shù)值,使得ITAE指標持續(xù)地滿足需求,即實際系統(tǒng)模型輸出與含待辨識參數(shù)的系統(tǒng)模型輸出不斷逼近一致。設定仿真時間截止或ITAE指標滿足要求時,可得出待辨識參數(shù)具體值。

3 模型參數(shù)辨識舉例驗證

3.1 實測數(shù)據(jù)獲取

由于電站實測數(shù)據(jù)較難獲取,因此本文用仿真方法獲取替代實測響應數(shù)據(jù)。采用Simulink建立仿真數(shù)據(jù)采集模塊。

用階躍響應外加一定量噪聲干擾當作水輪機輸入數(shù)據(jù),導入workspace2,獲取水輪機輸出響應數(shù)據(jù)(功率)替代實測輸出響應數(shù)據(jù)并導入workspace1。數(shù)據(jù)如圖3所示。

對理想水輪機模型進行辨識。理想水輪機數(shù)學模型如式(9)所示。在參數(shù)待辨識模型的辨識中采取的輸入輸出數(shù)據(jù)如圖3所示,對辨識后模型輸出和原真實系統(tǒng)模型輸出數(shù)據(jù)進行對比時,真實系統(tǒng)模型和辨識后系統(tǒng)模型均采用相同的階躍擾動,外附加小噪聲輸入激勵。

3.2 批處理最小二乘法參數(shù)辨識

運用最小二乘法,根據(jù)實測輸入輸出數(shù)據(jù)辨識Tw參數(shù)值。對理想水輪機傳遞函數(shù)模型進行拉普拉斯反變換,得到式(21)。

引入白噪聲序列 [xi],運用批處理最小二乘法,根據(jù)實測(樣本)數(shù)據(jù),對待辨識模型進行參數(shù)辨識并計算其誤差,對辨識結(jié)果進行評價。

輸入只給定階躍擾動的辨識后模型輸出響應與實測輸出響應關(guān)系,如圖4所示。

實測(樣本)與最小二乘辨識誤差曲線如圖5所示。

求得均方根誤差為0.034 4,平均絕對百分比誤差為0.076 9,且從圖5可以看出,辨識后響應和實測響應在一定程度上比較貼切,能在較短時間內(nèi)實現(xiàn)基本逼近。該方法不僅容易實現(xiàn)對線性系統(tǒng)模型參數(shù)的辨識,而且能達到較好的辨識效果。

3.3 基于OCD程序與ITAE準則的參數(shù)辨識

搭建可以實現(xiàn)參數(shù)辨識的Simulink模塊,如圖6所示。

其中將待辨識參數(shù)值設定為a1、a2。將模擬實際系統(tǒng)模型的輸出與含待辨識參數(shù)的模型輸出作為誤差給誤差積分準則ITAE模塊。

調(diào)用OCD程序,打開人機對話界面,如圖7所示。

輸入Simulink模塊名,輸入待辨識參數(shù)變量名,根據(jù)經(jīng)驗設定待辨識參數(shù)值上下界作為人為約束,設定仿真時長,點擊Create File按鈕,繼續(xù)連續(xù)點擊Optimize按鈕,即可讓Simulink模塊參與OCD程序參數(shù)辨識訓練,得到最終仿真結(jié)果和待辨識參數(shù)值。

得到a1、a2分別為-1.974 1、0.972 9,易分離得到Tw值為1.96。得到實際系統(tǒng)模型和辨識系統(tǒng)模型輸出曲線如圖8所示。

實際系統(tǒng)輸出與辨識系統(tǒng)輸出誤差曲線如圖9所示。

得到均方根誤差為0.008 3,平均絕對百分比誤差為0.018 4。可見誤差明顯低于最小二乘法辨識結(jié)果,辨識精確度更高。

將最小二乘法與基于OCD程序和ITAE準則的辨識結(jié)果進行比較,如圖10所示,從中可直觀地看出OCD辨識結(jié)果比最小二乘的辨識仿真結(jié)果更貼合實際(樣本值),但最小二乘辨識結(jié)果響應速度更快。

4 結(jié)語

本文探討了批處理最小二乘法與基于數(shù)值最優(yōu)化技術(shù)的最優(yōu)控制方法對實測數(shù)據(jù)(樣本)的參數(shù)辨識效果,通過參數(shù)辨識后得到的輸出功率與現(xiàn)場實測輸出(樣本)具有較高的擬合度。批處理最小二乘法能夠?qū)崿F(xiàn)較好的線性辨識效果,可較容易地獲取辨識后參數(shù)具體值,且操作簡便。基于數(shù)值最優(yōu)化技術(shù)的最優(yōu)控制方法辨識結(jié)果較最小二乘法具有更高的與實測數(shù)據(jù)擬合度。基于數(shù)值最優(yōu)化技術(shù)的最優(yōu)控制方法先對待辨識傳遞函數(shù)模型進行差分方程變換,而后運用批處理最小二乘法進行精細辨識,符合模型參數(shù)辨識數(shù)學步驟,再參考得到線性模型參數(shù)辨識的人機交互界面,避免底層程序編程,根據(jù)輸入輸出實測數(shù)據(jù)進行參數(shù)辨識。該方法具有一定的普適性和參考價值。

以上兩種辨識方法對基于實測(樣本)數(shù)據(jù)的線性系統(tǒng)模型參數(shù)辨識有良好的仿真結(jié)果,參數(shù)辨識值與原系統(tǒng)參數(shù)值逼近,驗證了其對給定實測輸入輸出數(shù)據(jù)的真實系統(tǒng)在小范圍波動下,在近似線性化模型參數(shù)精細化辨識上具有可行性。參數(shù)辨識后得到的模型可替代真機模型進行仿真研究工作,可得到符合實際系統(tǒng)期望的模型。下一步將對非線性系統(tǒng)模型參數(shù)辨識進行研究,驗證OCD人機界面應用于非線性系統(tǒng)辨識的可行性。

參考文獻:

[1] 陳志環(huán). 水輪機調(diào)節(jié)系統(tǒng)的參數(shù)辨識與控制策略研究[D]. 武漢:華中科技大學,2017.

[2] 丁聰,把多鐸,陳帝伊,等. 混流式水輪機調(diào)節(jié)系統(tǒng)的建模與非線性動力分析[J]. 武漢大學學報:工學版,2012,45(2):187-192.

[3] 彭伊麗. 用于水電控制仿真中的水輪機模型辨識[D]. 武漢:華中科技大學,2011.

[4] 華中科技大學. 一種水輪機調(diào)節(jié)系統(tǒng)的參數(shù)辨識方法:CN201510759863.1[P]. 2016-02-03.

[5] 高曉光,唐戢群,劉德發(fā), 等. 水輪機及調(diào)節(jié)系統(tǒng)模型辨識試驗研究[J]. 人民長江,2016,47(13):84-89.

[6] 方日升. 水輪機引水系統(tǒng)精細化模型參數(shù)辨識研究[J]. 電力與能源,2018,39(2):261-265.

[7] DEMELLO F P,KOESSLER R J,AGEE J,et al. Hydraulic turbine and turbine control models for system dynamic studies[J].? IEEE Transactions on Power Systems, 1992, 7(1):167-179.

[8] 張記坤,曾云,王芳芳,等. 混流式水輪機調(diào)節(jié)系統(tǒng)非線性建模與分叉分析[J]. 軟件導刊,2019,18(12):102-107.

[9] 張楚. 抽水蓄能機組調(diào)節(jié)系統(tǒng)非線性辨識與優(yōu)化控制研究[D]. 武漢:華中科技大學,2018.

[10] 吳凡,李偉雄. 基于MATLAB系統(tǒng)辨識工具的系統(tǒng)辨識[J]. 河北農(nóng)機,2016,(11):59-60.

[11] 陳嵐峰,張亞琴,程立英,等. 基于數(shù)據(jù)的MATLAB系統(tǒng)辨識工具箱模型識別[J]. 沈陽師范大學學報(自然科學版),2013,31(4):527-530.

[12] 龐中華,崔紅. 系統(tǒng)辨識與自適應控制MATLAB仿真(修訂版)[M]. 北京:北京航空航天大學出版社,2013.

[13] 薛定宇. 控制系統(tǒng)計算機輔助設計-MATLAB語言與應用(第3版)[M]. 北京:清華大學出版社,2012.

[14] 盧志強. 連續(xù)時間系統(tǒng)課程的離散化方法綜述[J]. 中國科教創(chuàng)新導刊,2010,(34):85-86.

[15] 王貞祥,高立群,李瑋,等. 系統(tǒng)辨識與參數(shù)估計[M]. 沈陽:東北大學出版社,1993.

[16] 李偉,王洪民,唐崢. 基于遞推最小二乘法的轉(zhuǎn)向系統(tǒng)參數(shù)辨識[J]. 重慶交通大學學報(自然科學版),2019,38(8):124-128.

[17] 張國洋. Matlab在水輪機模型辨識中的應用[J]. 自動化應用,2016(1):19-20,26.

[18] 黃莉,李咸善. 基于Simulink的水輪機動態(tài)仿真模型研究[J]. 三峽大學學報(自然科學版),2007,29(1):25-28.

[19] 吳卓璠. 水輪機調(diào)節(jié)系統(tǒng)非線性建模與參數(shù)辨識[D]. 武漢:華中科技大學,2016.

[20] 劉昌玉,梁學磊. 水輪機調(diào)節(jié)系統(tǒng)被控對象模型辨識[J]. 水電能源科學,2007,25(2):77-79.

[21] 柳海生,曾江,黃海穎,等. 抽水蓄能電站水泵水輪機效率計算方法[J]. 廣東電力,2014(9):1-5.

(責任編輯:江 艷)

主站蜘蛛池模板: 日韩欧美国产精品| 亚洲精品自在线拍| 欧美日本不卡| 久久这里只有精品国产99| 国产网站免费观看| 欧美、日韩、国产综合一区| 成人免费网站久久久| 国产在线精品香蕉麻豆| 精品国产免费观看| 国产成人艳妇AA视频在线| 国产探花在线视频| 日韩免费中文字幕| 亚洲一区二区在线无码| 国产激情第一页| 国产菊爆视频在线观看| 99re66精品视频在线观看| 精品国产网站| 欧美亚洲国产视频| 久久伊人色| 亚洲免费福利视频| 国产永久在线视频| 欲色天天综合网| 中国一级毛片免费观看| 在线五月婷婷| 国产精品久久久久久影院| 亚洲一区波多野结衣二区三区| 亚洲狠狠婷婷综合久久久久| 人人91人人澡人人妻人人爽| 欧美在线伊人| 91在线高清视频| 午夜毛片福利| 大香网伊人久久综合网2020| 久久久久久午夜精品| 亚洲欧美在线看片AI| 亚洲IV视频免费在线光看| 制服丝袜 91视频| 日韩小视频在线观看| 欧美一道本| 国产精品分类视频分类一区| 亚洲色图综合在线| 国产情侣一区| 色吊丝av中文字幕| 国产在线日本| 国产日韩精品一区在线不卡| 国产精品第一区在线观看| 高清乱码精品福利在线视频| 亚洲91精品视频| 亚洲精品自产拍在线观看APP| 国产精品lululu在线观看| 麻豆国产精品视频| 国产亚洲欧美在线人成aaaa| 激情成人综合网| 国产资源免费观看| 波多野结衣久久高清免费| 亚洲精品国产日韩无码AV永久免费网| 中文字幕2区| 直接黄91麻豆网站| 欧美日韩国产一级| 伊人色天堂| 国产女人18毛片水真多1| 午夜日b视频| 成人免费黄色小视频| 中文字幕在线视频免费| 久久青草免费91观看| 国产永久在线视频| 26uuu国产精品视频| 久久综合激情网| 国产精品妖精视频| 亚洲狠狠婷婷综合久久久久| 国产区在线观看视频| 97se亚洲综合在线韩国专区福利| 亚洲乱亚洲乱妇24p| 91九色国产在线| 六月婷婷精品视频在线观看 | www中文字幕在线观看| 自拍欧美亚洲| 亚洲一道AV无码午夜福利| 久久国产乱子伦视频无卡顿| 国产成人啪视频一区二区三区 | 亚洲欧美日韩中文字幕一区二区三区 | 国产免费羞羞视频| 伊人久久婷婷|