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

基于可變設計參數的平面Acrobot位置快速控制方法

2017-09-22 01:18:41王亞午賴旭芝吳敏
電機與控制學報 2017年9期
關鍵詞:系統設計

王亞午,賴旭芝,吳敏

(1.中國地質大學(武漢) 自動化學院,湖北 武漢 430074;2.復雜系統先進控制與智能自動化湖北省重點實驗室,湖北 武漢 430074)

基于可變設計參數的平面Acrobot位置快速控制方法

王亞午1,2,賴旭芝1,2,吳敏1,2

(1.中國地質大學(武漢) 自動化學院,湖北 武漢 430074;2.復雜系統先進控制與智能自動化湖北省重點實驗室,湖北 武漢 430074)

為平面Acrobot系統(無重力項)提出一種基于可變設計參數的快速位置控制方法。首先,根據平面Acrobot系統的完全可積特性,獲取驅動關節與欠驅動關節的角速度和角度約束關系。其次,根據系統的物理結構,推導出末端點坐標位置和兩關節角度之間另一類約束關系。然后,基于這兩類約束關系,利用粒子群優化算法求取目標位置對應的兩桿目標角度。接著,通過構造李雅普諾夫函數為平面Acrobot系統設計控制器,并根據控制器中設計參數的特性,選擇一個可變的設計參數,減少控制系統到達目標的時間。利用連桿角速度和角度約束關系,通過控制驅動桿到目標角度同時連帶控制欠驅動桿到目標角度,實現系統的控制目標。最后,仿真結果驗證所設計控制方法的有效性和快速性。

平面欠驅動機械系統;平面Acrobot;位置控制;粒子群優化算法;非線性函數

0 引 言

欠驅動機械系統是指控制輸入個數少于系統自由度個數的非線性系統[1-3],它能夠利用較少的驅動裝置完成控制任務,具有重量輕、低能耗、靈活等優點;然而,驅動裝置的減小使得該類系統控制器設計的難度增加,很多適用于全驅動系統的控制方法無法直接推廣到這類系統中。而目前,國內外對垂直欠驅動機械系統(含重力項)的研究較為成熟[4-6],但對平面欠驅動機械系統(不含重力項)的研究還在起步階段。

針對垂直欠驅動機械系統,現在普遍采用的控制方法是分區控制方法[1,7-8],即將控制空間劃分為搖起區和平衡區。針對搖起區的控制問題,主要采用基于系統模型[9]、運動姿態[10]、系統能量[11]和魯棒性能[12]等指標的智能控制方法;而平衡區的控制主要通過在平衡點附近將系統近似線性化,并采用線性二次調節器(LQR)來實現系統的穩定控制[13]。而對于平面欠驅動機械系統的控制問題,由于其動力學模型不含重力項,導致其在水平面內的任意靜止姿態都為它的平衡點,并且其在平衡點附近的線性近似模型不滿足可控條件[14],無法采用垂直欠驅動系統中已成熟的LQR方法實現系統穩定控制。

早期研究中學者們為降低平面欠驅動系統的控制難度,在模型中考慮摩擦力[15-16]或者在欠驅動關節加裝制動裝置或彈簧[17-18],以此完成平面欠驅動系統的控制任務;但這些方法改變了系統的欠驅動性質。文獻[19]通過為平面Pendubot離線規劃一條軌跡,并利用滑模控制技術實現系統軌跡的穩定跟蹤。文獻[20]針對平面Pendubot的控制問題,通過分解系統的運動,利用分層模糊控制器實現系統的控制。但這些研究大多針對第一關節為驅動的平面機械系統,而對于第一關節為欠驅動的平面機械系統的研究則較少。Oriolo等人[21]對平面欠驅動系統的可積性進行了初步研究,并證明了第一關節為欠驅動的平面Acrobot具有完全可積分的特性。之后,文獻[22]基于平面Acrobot的完全可積特性,提出了基于軌跡特性的平面Acrobot運動控制方法,但從其仿真結果來看,系統各個狀態到達穩定的時間較長,并且沒有提出基于末端點位置獲取連桿角度的方法。

基于文獻[21]和文獻[22]的研究,為了克服系統控制時間較長的問題,提出一種快速實現末端點位置控制的方法。首先,基于平面Acrobot的完全可積特性,獲取驅動桿和欠驅動桿的角速度約束關系和角度約束關系,即通過控制驅動桿(第二連桿)實現欠驅動桿(第一連桿)的連帶控制。然后,根據系統的物理結構,獲取末端點位置和兩桿角度之間的約束關系。基于系統這兩類約束關系,利用粒子群優化(PSO)算法求取末端點目標坐標位置對應的兩桿目標角度值。最后,通過構造李雅普諾夫函數設計控制器,并根據控制器中設計參數的特性,將設計參數設計成一類非線性函數,減少控制系統到達目標的時間。通過控制驅動桿達到目標角度同時連帶控制欠驅動桿達到目標角度,實現系統末端點位置控制任務。仿真結果驗證所提控制方法不僅能夠完成系統末端點的位置控制任務,并且控制時間較短。

1 系統模型

第一關節為欠驅動的平面兩連桿系統(平面Acrobot)的物理結構如圖1所示:

圖1 平面Acrobot系統Fig.1 Planar Acrobot system

對于第i連桿(i=1,2):mi表示質量;Li表示長度;li表示第i連桿的質心到第i關節的長度;Ji表示轉動慣量;q1為第1連桿相對于坐標y軸方向的角度,q2為第2連桿相對于第1連桿延長線方向的角度;(x,y)為系統末端點坐標位置。

根據歐拉-拉格朗日公式可得平面Acrobot系統的動力學方程為

(1)

上述各元素的具體形式如下:

M11(q)=a1+a2+2a3cosq2,

M12(q)=M21=a2+a3cosq2,

M22(q)=a2,

其中a1、a2、a3為系統結構參數,具體表達式如下:

a3=m2L1l2。

由于平面Acrobot系統運動時處于水平面內,不含重力項,因而簡化了系統的動力學模型;但重力項的缺失會使得系統在水平面內的任意靜止姿態都成為它的平衡點,并且也導致在平衡點附近的線性近似模型不滿足線性可控條件[15],故無法采用垂直欠驅動系統中已成熟的LQR方法實現系統的穩定控制。因此,為實現平面Acrobot系統從任意初始位置到目標位置的控制,下面將基于平面Acrobot的完全可積特性獲取系統驅動桿與欠驅動桿之間的狀態約束關系,并基于此約束關系提出平面Acrobot的位置控制方法。

2 約束關系分析

由動力學方程(1)可知,系統的二階微分約束為

(2)

系統慣性矩陣的欠驅動部分Mu(q)=[M11(q),M12(q)]對時間求導,可得

(3)

(4)

將式(4)代入到式(2)中得

(5)

式(5)對時間積分得

(6)

(7)

從式(7)可知,當第二桿的角速度為0時,第1桿的角度速度也為0,即只要能控制第2桿到靜止狀態,平面Acrobot系統就能穩定。

而M11(q)是正定對稱矩陣M(q)對角線上的元素,則有M11(q)>0。又因M(q)不包含q1,所以式(7)可寫成

(8)

式(8)對時間積分,可得

(9)

其中A1和A2是與系統結構參數有關的量,具體形式分別為:

(10)

其中floor(·)為向下取整函數。

1)k≥0時,

f(q2)-f(q2(0))+kπ

f(q2)-f(q2(0))+kπ

(11)

其中,c=-A1f(q2(0))+q1(0)+q2(0)/2。

從角度約束關系式(11)可知,當控制驅動桿(第2桿)到達某個角度,欠驅動桿(第1桿)也將被連帶控制到某個角度,即可通過控制驅動桿來實現欠驅動桿的連帶控制。下面,通過分析系統的幾何結構,同時結合角度約束關系式(11),來獲取系統末端點目標坐標位置對應的兩桿目標角度。

3 目標角度求解

在討論目標角度求解之前,首先需建立末端點的坐標位置(x,y)和兩桿角度q1,q2之間的關系,如圖2所示。

圖2 平面Acrobot物理結構簡圖Fig.2 Physical structure diagram of planar Acrobot

根據平面Acrobot的物理結構簡圖,運用幾何知識可獲得末端點坐標位置與兩桿角度之間的約束關系為:

(12)

因此,當系統末端點的目標坐標位置設計為(xd,yd)時,為尋找同時滿足約束關系式(11)和式(12)的目標角度q1d和q2d,首先將式(11)代入式(12)中得:

(13)

然后根據式(13)可定義如下目標函數:

h(q2)=|x-xd|+|y-yd|。

(14)

可知,h(q2)=0的解就是驅動桿目標角度q2d。當將q2d代入約束關系式(11)中,即可求得對應的欠驅動桿目標角度q1d。

考慮到直接求解h(q2)=0較為困難,并且在誤差允許的范圍內(h(q2)≤e1),目標坐標位置對應的目標角度具有多解性。因此,這里選取粒子群優化算法作為求解方法。記sj=q2為粒子的位置。從而,第j(j=1,2,…,n)個粒子的進化方程為:

(15)

其中:w為慣性權重;c1和c2為權重因子;r1和r2為在[0,1]當中的隨機值;sj(K)表示第j個粒子的位置;vj(K)表示第j個粒子的速度;pj代表第j個粒子的個體歷史最優位置;pg表示群體歷史最優位置;K為粒子進化代數;n為粒子群的大小。

個體歷史最優位置更新規則:

(16)

群體歷史最優值更新規則

pg=argmin{h(pj)}。

(17)

為防止粒子的進化速度過大,這里加入速度限制方程,如下式所示:

|vj(K+1)|≤vmax。

(18)

其中vmax為正的常數。

求解兩桿目標角度的算法步驟如下:

1)初始化n個粒子的初始位置sj(0)與初始速度vj(0)及個體歷史最優位置pj與群體歷史最優位置pg。

2)按照粒子的進化方程式(15)來更新各個粒子的位置和速度,由式(18)來限制粒子進化速度。當進化代數K=Kmax,退出程序。

3)將新一代n個粒子的位置代入到目標函數h(·)中,求得此新一代n個粒子的適應度值。

4)根據各個粒子的適應度值,按照式(16)和式(17)來更新個體歷史最優位置pj和群體歷史最優位置pg。

5)當粒子群中最大的適應度值小于或等于e1時,群體歷史位置pg即為驅動桿目標角度q2d,程序轉到第六步;否則,轉到第2步。

6)將q2d代入到公式(11)中,求出相應欠驅動桿目標角度q1d,退出程序。

當控制驅動桿到達目標角度q2d同時連帶欠驅動桿到目標角度q1d時,系統末端點運動到目標坐標位置(xd,yd)。基于此思想,下面設計控制器實現平面Acrobot末端點的位置控制。

4 控制器設計

從上面分析可知,當控制驅動桿到達目標角度q2d時,欠驅動桿也將被連帶控制達到目標角度q1d,從而實現系統的控制目標。本節將依據系統這一特性,通過為平面Acrobot系統構造李雅普諾夫函數來設計控制器,并針對控制器的設計參數,提出改進計劃,實現系統末端點位置的快速控制。

(19)

其中,Fi(X)、Gi(X)是關于X的非線性函數,具體形式如下:

構造如下李雅普諾夫函數為

(20)

其中,x2d=q2d。

式(20)對時間求導得

(21)

(22)

若設計

x2-x2d+F2(X)+G2(X)τ2=-rx4。

(23)

其中,r(r>0)為設計參數。

則將式(23)代入式(22)中有

(24)

因此,控制器可設計為

(25)

(26)

將式(25)或式(26)代入式(19)中可得到以下閉環系統為

(27)

閉環系統式(27)的框圖如圖3所示。

圖3 閉環系統(27)的框圖Fig.3 Block diagram of closed system(27)

Φ={X∈R4|0≤V(X)≤C}。

(28)

其中C為正的常數。

系統式(27)任意始于Φ內的解X對于所有t≥0仍處于Φ。記Ψ為系統(27)的不變集,

(29)

(30)

(31)

由LaSalle不變原理[23]可知,當t→,系統式(27)始于Φ內的解X趨于M,即:

因此,系統式(19)在控制器式(25)作用下可以運動并穩定到目標狀態,即平面Acrobot系統末端點到達所設定的目標坐標位置(xd,yd)。

從圖3可以看出,所設計控制器中包含1個PD控制器,其比例系數為1,微分時間參數為設計參數r(r>0)。當r較小時,e變化快(即x2變化快),超調大;而當r較大時,e變化慢(即x2變化慢),超調小。基于此想法,當第2桿角度x2與目標角度x2d偏差絕對值|e|較大時,調小r值,使e變化快;當|e|較小時,調大r值,防止系統超調。而且,為不影響系統的收斂性,應保證r>0。

基于上述思想,可將r設計成與驅動桿角度相關的一類非線性函數,即|x2d-x2|較大時,r取小;|x2d-x2|較小時,r取大。這里以分段函數和高斯函數為例進行說明。

1)r為分段函數,如下式所示:

(32)

其中:ed、rmax和rmin(rmax>rmin)為正的常數。

式(32)的函數圖形如圖4所示。

圖4 分段函數Fig.4 Piecewise function

由圖4可知,當x2離目標角度x2d較遠時,r=rmin;當x2離目標角度x2d較近時,r=rmax。并且r一直處于正值,不影響系統式(27)的穩定性。那么控制器式(25)變為:

(33)

2)r為高斯函數,如下式所示:

(34)

其中A和σ都是正的常數。

式(34)的函數圖形如下

圖5 高斯函數Fig.5 Gauss function

由圖5可知,當第2桿角度x2離目標角度x2d較遠時,采用較小的r;而隨著第2桿角度x2接近目標角度x2d時,r也隨之增大。并且r一直處于正值,不影響系統式(27)的穩定性。因此控制器式(25)變為

(35)

5 仿真結果

按照上述設計方案,利用Matlab環境下的SIMULINK平臺搭建系統仿真模型,完成平面Acrobot系統仿真控制。平面Acrobot系統的仿真參數如表1所示。

表1 平面Acrobot系統的仿真參數Table 1 Simulation parameter of plannar Acrobot system

圖6 r為常數時的仿真結果(r=1.8)Fig.6 Simulation results when r is a constant

1)r為常數,仿真結果如圖6所示。

從圖6(a)可看出,兩桿均達到目標角度;從圖6(c)可看出,系統末端點在6 s時已經達到所設計的目標坐標位置,說明所設計的控制器有效。

2)r為分段函數,仿真結果如圖7所示。

從圖7(a)可看出,兩桿均達到目標角度;從圖7(c)可看出系統末端點在2.7 s時已經達到目標位置(1.25,1.55)。說明所設計的控制器不僅有效,而且可縮短控制系統達到目標的時間。

3)r為高斯函數,仿真結果如圖8所示。

圖8 r為高斯函數時的仿真結果(A=60,σ=1)Fig.8 Simulation results when r is a gauss function

從圖8(a)可看出,兩桿均達到目標角度;從圖8(c)可看出系統末端點在2 s時已經達到目標位置(1.25,1.55)。同樣說明,所設計的控制器不僅有效,并且可縮短控制系統達到目標的時間。

因此,所設計的控制器不但能夠完成系統末端點的位置控制任務,而且當將控制器的設計參數r設計成分段函數或者高斯函數這一類非線性函數時,可以較大地減少控制系統末端點到達目標坐標位置的時間。

6 結 論

為實現平面Acrobot系統末端點從初始位置到目標位置的快速控制,本文提出了一種基于可變設計參數的控制方法。該方法利用平面Acrobot系統的完全可積性質,通過控制驅動桿實現了欠驅動桿的連帶控制。并根據平面Acrobot系統的兩類約束關系定義目標函數,由PSO算法解出兩桿目標角度值。最后通過構造李雅普諾夫函數來設計控制器,實現了平面Acrobot末端點的位置控制任務。特別是根據控制器中設計參數r的特性,將r設計成與驅動桿角度相關的一類非線性函數,實現了系統末端點位置的快速控制。仿真結果的對比驗證了所提控制方法可以較大地減少系統到達控制目標的時間。

[1] LAI X Z,PAN C Z,WU M,et al.Unified control of n-link underactuated manipulator with single passive joint:A reduced order approach[J].Mechanism and Machine Theory,2012,56(12):170.

[2] 張靜,孫強,關婷婷.TORA系統的自適應神經網絡輸出反饋控制[J].哈爾濱理工大學學報,2015,20(6):83.ZHANG Jing,SUN Qiang,GUAN Tingting.Output feedback control of the TORA system based on adaptive neural network[J].Journal of Harbin University of Science and Technology,2015,20(6):83.

[3] 高丙團,陳宏鈞,張曉華.欠驅動機械系統控制設計綜述[J].電機與控制學報,2006,10(5):541.GAO Bingtuan,CHEN Hongjun,ZHANG Xiaohua.Control design for underactuated mechanical systems:a survey[J].Electric Machines and Control,2006,10(5):541.

[4] 劉殿通,易建強,譚民.不確定欠驅動非線性系統的模糊滑模控制[J].電機與控制學報,2003,7(3):215.LIU Diantong,YI Jianqiang,TAN Min.Fuzzy sliding mode control of uncertain underactuated nonlinear systems[J].Electric Machines and Control,2003,7(3):215.

[5] LAI X Z,ZHANG A C,SHE J H,et al.Motion control of underactuated three-link gymnast robot based on combination of energy and posture[J].IET Control Theory and Application,2011,5(13):1484.

[6] 武俊峰,王勝達.模糊控制在二級倒立擺系統中的應用[J].哈爾濱理工大學學報,2013,18(4):12.WU Junfeng,WANG Shengda.The application of fuzzy control in double inverted pendulum[J].Journal of Harbin University of Science and Technology,2013,18(4):12.

[7] XIN X,SHE J H,YAMASAKI T,et al.Swing-up control based on virtual composite links for n-link underactuated robot with passive first joint[J].Automatica,2009,45(9):1986.

[8] ZHANG M J,TARNTJ.Hybrid control of the pendubot[J].IEEE/ASME Transactions on Mechatronics,2002,7(1):79.

[9] PATRICIO O J,OMAR S J,VIRGILIO L M.Toward a generalized sub-optimal control method of underactuated systems[J].Optimal Control Applications & Methods,2012,33(3):338.

[10] XIN X,MITAT,KANEDAM.The posture control of a two-link free flying acrobot with initial angular momentum[J].IEEE Transactions on Automatic Control,2004,49(7):1201.

[11] FANTONI I,LOZANOR,SPONG M W.Energy based control of the pendubot[J].IEEE Transactions on Automatic Control,2002,45(4):725.

[12] LAI X Z,PAN C Z,WU M,et al.Robust stabilization and disturbance attenuation for a class of underactuated mechanical systems[J].Journal of Central South University,2012,19(9):2488.

[13] MARTON L,HODEL A S,LANTOS B,et al.Underactuated robot control:comparing LQR,subspace stabilization,and combined error metric approaches[J].IEEE Transactions on Industrial Electronics,2008,55(10):3724.

[14] LUO Yibiao,LAI Xuzhi,WU Min.A positioning control strategy of planar Acrobot[C].Proceeding of the 30th Chinese Control Conference,Yantai,China,2011:743-747.

[15] SUZUKI T,MIYOSHI W,NAKAMURA Y.Control of 2R underactuated manipulator with friction[C].Proceedings of the 37th IEEE Conference on Decision and Control,Tampa,USA,1998,2007-2012.

[16] MAHINDRAKAR A D,SHODHAN R,BANAVAR R N.Point-to-point control of a 2R planar horizontalunderactuated manipulator[J].Mechanism and Machine Theory,2006,41(7):838.

[17] ARAIH,TACHIS.Position control of manipulator with passive joints using dynamic coupling[J].IEEE Transactions onRobotics and Automation,1991,7(4):528.

[18] POPESCU C,WANG Y,ROTH Z.Passivity based control of spring coupled underactuated horizontaldouble pendulum[C].Proceedings of the Florida Conference of recent Advances in Robotics,Florida,USA,2003:1-17.

[19] ALFREDO J R F,JAIME A G,Rafael C L.Trajectory planning and control of an underactuated planar2R manipulator[C].Proceedings of the IEEE International Conference on Control Applications,Mexico,2001:548-552.

[20] 劉慶波,余躍慶.基于遺傳算法的欠驅動機器人模糊控制器設計[J].系統仿真學報,2008,20(8):2097.LIU Qingbo,YU Yuejin.Fuzzy controller design for underactuated robots based on GA method[J].Journal of System Simulation,2008,20(8):2097.

[21] ORIOLOG,NAKAMURAY.Control of mechanical systems with second-order nonholonomic constraints:Underactuated manipulators[C].Proceedings of 30th IEEE Conference on Decision and Control,Brighton,England,1991:2398-2403.

[22] CAO Shengqiang,LAI Xuzhi,WU Min.Motion control method of planar Acrobot based on trajectory characteristics[C].Proceeding of the 31th Chinese Control Conference,Hefei,China,2012:4910-4915.

[23] KHALIL H K.Nonlinear systems (3rd Ed.)[M].New Jersey:Prentice Hall,2002.

(編輯:張 楠)

RapidpositioncontrolapproachbasedonvariabledesignparameterforplanarAcrobot

WANG Ya-wu1,2,LAI Xu-zhi1,2,WU Min1,2

(1.School of Automation,China University of Geosciences,Wuhan 430074,China;2.Hubei key Laboratory of Advanced Control and Intelligent Automation for Complex Systems,Wuhan 430074,China)

A rapid position control approach based on a variable design parameter was presented for a planar Acrobot which has no gravity item.Firstly,angle constraint relationship and angular velocity constraint relationship between active joint and passive joint were obtained by employing the complete integral characteristics of the planar Acrobot.Next,another kind of the constraint relationship between the endpoint position and two angles of joint was deduced according to the physical structure of the planar Acrobot.Then,the target angles of joint were obtained by using particle swarm optimization (PSO) algorithm based on the above two kinds of constraint relationship.Finally,a controller was designed for planar Acrobot by constructing a Lyapunov function,and a variable design parameter of the controller was chosen to reduce the control time according to the characteristics of parameter.When the active link was controlled to the target angle,the passive link was also controlled jointly to its target angle by utilizing angular velocity constraint relationship and angle constraint relationship.That is,the endpoint of the planar Acrobot was driven to the target position.Simulation results demonstrate the effectiveness and rapidity of the proposed control approach.

planar underactuated mechanical system; planar Acrobot; position control; particle swarm optimization algorithm; nonlinear function

10.15938/j.emc.2017.09.015

TP 4

:A

:1007-449X(2017)09-0110-09

2015-01-20

國家自然科學基金(61374106);湖北省自然科學基金(2015CFA010);“111計劃”項目(B17040)

王亞午(1990—),男,博士研究生,研究方向為欠驅動機器人控制、非線性系統控制;賴旭芝(1966—),女,教授,博士生導師,研究方向為智能控制、機器人控制和非線性系統控制;吳 敏(1963—),男,教授,博士生導師,研究方向為魯棒控制、智能控制和過程控制。

賴旭芝

猜你喜歡
系統設計
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
何為設計的守護之道?
現代裝飾(2020年7期)2020-07-27 01:27:42
《豐收的喜悅展示設計》
流行色(2020年1期)2020-04-28 11:16:38
基于PowerPC+FPGA顯示系統
半沸制皂系統(下)
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
主站蜘蛛池模板: 欧美亚洲日韩不卡在线在线观看| 性喷潮久久久久久久久| 亚洲黄色成人| 亚洲国产av无码综合原创国产| 午夜性刺激在线观看免费| 亚洲男人在线| 亚洲精品福利视频| 97亚洲色综久久精品| 国产一区自拍视频| 人妻精品久久无码区| 日韩av无码精品专区| 国产免费精彩视频| 国产成人高清在线精品| 亚洲中久无码永久在线观看软件| 亚洲综合片| 在线无码av一区二区三区| 精品久久久久久中文字幕女| 午夜视频免费试看| 国产欧美日韩精品第二区| 国产精品亚洲日韩AⅤ在线观看| 成人福利在线看| 亚洲一区国色天香| 久久99国产综合精品1| 亚洲欧美不卡中文字幕| 激情亚洲天堂| 国产午夜人做人免费视频| 丝袜久久剧情精品国产| 日本高清免费不卡视频| 五月婷婷亚洲综合| 国产97公开成人免费视频| 国产亚洲精品在天天在线麻豆 | 日本免费精品| 亚洲精品卡2卡3卡4卡5卡区| 亚洲国产精品日韩欧美一区| 亚洲国产午夜精华无码福利| 久久99国产综合精品女同| 成人精品午夜福利在线播放| 日本免费一区视频| 97人人做人人爽香蕉精品| 成人年鲁鲁在线观看视频| 久久无码免费束人妻| 性激烈欧美三级在线播放| 亚洲视频免费在线看| 亚洲人在线| 亚洲欧洲日产国产无码AV| 欧美另类精品一区二区三区| 国产成人无码AV在线播放动漫 | 亚洲视屏在线观看| 波多野结衣的av一区二区三区| 国产成a人片在线播放| 亚洲大尺度在线| 欧美精品导航| 亚洲欧美日本国产综合在线| 国产成人免费| 亚洲精品无码av中文字幕| 国产成人精品亚洲日本对白优播| 国内丰满少妇猛烈精品播 | 国产精品漂亮美女在线观看| 国产一级二级在线观看| 中文字幕有乳无码| 亚洲第一黄色网| 午夜一级做a爰片久久毛片| 狠狠综合久久| 熟妇无码人妻| 欧美成人h精品网站| 国产后式a一视频| 久久99精品久久久大学生| 波多野结衣第一页| 毛片网站观看| 国产免费久久精品99re丫丫一| 欧美国产综合视频| 69国产精品视频免费| 一本久道热中字伊人| 国产欧美亚洲精品第3页在线| 亚洲永久免费网站| 亚洲va精品中文字幕| 国产精品午夜福利麻豆| 国产福利在线免费观看| 99re在线视频观看| 亚洲第一色视频| av一区二区无码在线| 国产91在线|中文|