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

考慮預測響應值波動的多響應優化設計

2018-07-27 03:09:36汪建均屠雅楠
系統工程與電子技術 2018年8期
關鍵詞:優化方法模型

汪建均, 屠雅楠

(南京理工大學經濟管理學院, 江蘇 南京 210094)

0 引 言

近年來,隨著產品設計的復雜化以及顧客需求層次的多樣化,在產品或工藝過程的優化設計中往往需要同時考慮多個質量特性,因此多響應優化設計在持續性質量改進活動中顯示出越來越重要的地位和作用[1]。多響應優化設計旨在尋求一組最佳的參數設計值,使多個質量特性能夠同時達到最優的設計值。然而,在多響應優化設計的建模過程中,往往無法準確估計模型的回歸系數。如果忽略模型參數的這種不確定性,試驗結果將會產生誤差,甚至得出不科學的結果[2]。此外,在產品的生產過程中,存在多種噪聲因子的干擾,這些噪聲因子會對質量特性的輸出結果造成很大影響,使得預測響應值存在較大的波動,導致試驗結果難以有效地復現[3]。

目前解決多響應優化設計問題的主要方法有以下幾種:滿意度函數法、馬氏距離法、多元質量損失函數方法、貝葉斯后驗概率法等[4]。滿意度函數法最早是由文獻[5]在1965年提出的,經過眾多質量學者的不斷改進和完善,形成了一系列改進的滿意度函數。其中,最常用的一種是由文獻[6]所提出的綜合滿意度函數。該滿意度函數法操作簡單、易于理解,是使用最廣泛的多響應優化設計方法。然而,傳統的滿意度函數往往忽略了響應之間的相關性以及噪聲因子對優化結果的影響[7],對響應的預測效果以及優化結果的可靠性考慮較少[8-9]。在這種情況下,通過上述滿意度函數法所得到的優化結果往往是不可靠的,甚至會產生錯誤的結果。文獻[10]提出了一種指數滿意度函數,利用擬合優度提高了優化結果的可靠性,同時通過調節滿意度函數的形狀考慮了模型的預測能力。但是,該方法沒有考慮響應之間的相關性。文獻[11]提出一種改進的雙指數滿意度函數,該方法引入了質量損失函數來考慮響應間的相關性。但是,該方法也存在一些不足,未考察響應的預測誤差,同時響應的質量損失難以有效確定。文獻[12]提出一種改進的滿意度函數法,其主要優勢在于將雙響應曲面方法引入到滿意度函數中,能夠有效解決均值響應與方差響應之間的相關性問題。文獻[7]提出了改進的滿意度函數模型,該方法同時將可控因子和噪聲因子引入到滿意度函數,考慮了噪聲因子對多響應優化結果的影響。馬氏距離方法作為多響應優化設計的一種常用方法,也引起廣大質量學者的廣泛關注和重視。例如,文獻[13]提出的馬氏距離只考慮了響應的方差-協方差結構,未考察各個響應的相對重要性。為此,文獻[14]提出了加權馬氏距離,使優化結果具有一定的可靠性。該方法既考慮了響應的方差-協方差結構,還考慮了各個響應在經濟上的相對重要性。但是,該方法中的相對重要性受主觀因素的影響。如果權重設置不合理,該方法的優化效果會受到較大的影響。此外,馬氏距離方法主要是通過響應之間的方差-協方差結構來考慮多響應之間的相關性,但是通常忽視了響應的預測能力和優化結果的可靠性。多元質量損失函數也是目前最為常用的多響應優化設計方法之一,也引起眾多質量研究者的廣泛興趣。例如,文獻[15]提出的多元損失函數,考慮了過程的經濟性和響應間的相關性。然而,該方法存在較大的局限性,僅僅適合解決望目質量特性的問題。為此,文獻[16]擴展了文獻[15]的研究工作,將該方法的適用范圍擴大到望大質量特性和望小質量特性問題。文獻[17]提出的多元質量損失函數,考慮了響應的預測能力,但是忽略了響應之間的相關性。因此,文獻[18]結合文獻[15]和文獻[17]提出的兩種方法,構建了新的多元質量損失函數,該方法能夠同時考慮響應之間的相關性以及響應的預測能力。但是,該方法仍然存在缺陷,最主要的問題是無法評估優化結果的可靠性。文獻[19]提出一種新的質量損失函數,該方法能夠合理測算響應偏離目標值的程度。文獻[20]提出一種加權的質量損失函數,該方法同時考慮了二次損失的期望和方差。文獻[21]提出一種動態多響應穩健優化模型,該模型將多元偏正態分布和響應曲面法結合起來,建立了基于多元偏正態分布的期望損失函數,最后利用遺傳算法進行優化求解。隨著生產工藝復雜程度的日益提高,整個生產過程存在各種不確定性因素,盡管上述模型的理論基礎很完善,但是在各自的建模過程中沒有考慮模型參數不確定性和噪聲因子的共同影響,這些因素將使得預測響應值產生較大的波動,因此上述模型難以精確地描述整個生產過程。

針對上述問題,一些學者提出了各自的研究方法。文獻[22]提出,如果在建模過程中忽視模型參數的不確定性,該模型可能無法得到合理的優化結果。文獻[23]提出了一種貝葉斯后驗概率方法,該方法既考慮了過程分布的變化、試驗數據間的相關性,同時還考慮了模型參數的不確定性,然后通過蒙特卡羅模擬方法計算多個響應的后驗抽樣值落在規格限內的概率。文獻[24]考察了模型參數存在不確定性的多響應優化問題,結合貝葉斯方法開展研究工作。文獻[25]將帕累托最優策略引入多響應優化設計,該方法考慮了模型參數的不確定性問題,其主要特點在于能夠有效地利用圖形展示這種不確定性對多響應優化結果和用戶決策的影響。文獻[9]提出利用置信區間優化多響應問題,當模型參數存在不確定性時,該方法能夠給出更為合理的優化結果。文獻[2]提出一種基于置信區間的多元質量特性滿意參數設計方法,該方法考慮了模型參數不確定性的影響,利用雙響應曲面模型和置信區間求解滿足約束的相容性解集。文獻[26]提出一種新的損失函數方法,該方法不僅考慮了模型參數的不確定性,還引入了試驗過程中存在的實施誤差,并通過實例證明該方法能夠提高優化結果的穩健性和可靠性。另外,實際生產過程中存在各種噪聲因子,這些噪聲因子使得實際的生產過程難以通過單一的響應模型進行精確的刻畫。正如著名的統計學家Box所說的那樣“所有的模型都是錯誤,但是有些模型是有用的”。因此,一些研究人員開始關注如何通過組合建模方法以減少模型不確定性對優化結果的影響。例如,文獻[27]提出一種基于模型不確定性的響應曲面建模方法,該方法將包容性檢驗引入組合建模方法中,以解決模型不確定性問題。文獻[28]構建了一種新的損失函數,該方法結合位置和散度參數進行建模,并考慮了模型不確定性對試驗結果的影響。文獻[29]利用貝葉斯先驗信息,評估預測響應的后驗分布,該方法考慮了建模過程中可能存在的各種噪聲因子對優化結果的影響。文獻[30]通過置信區間衡量預測響應波動對優化結果的影響,該方法從置信區間端點選取最優值作為響應預測值,并利用圖形顯示預測響應波動的情況下優化結果的變化。文獻[31]提出一種基于貝葉斯方法的遞階優化算法,該方法利用滿意度函數求解初始解,再根據貝葉斯方法分析解的穩健性,綜合考慮了過程波動(可控因子波動或存在噪聲因子)和各種誤差的影響。文獻[32]結合質量損失函數和貝葉斯后驗概率方法提出了一種新的多響應優化方法。該方法在多響應優化設計中同時考慮多響應之間的相關性、模型參數的不確定性以及優化結果的可靠性。然而,上述方法未能進一步考慮噪聲因素以及模型參數的不確定性所導致的響應預測值波動對優化結果的影響,即不同研究方法所獲得的優化結果在試驗可重復性上的差異。在多響應優化設計中,如果忽略模型參數的不確定性,將難以精確地獲得生產過程輸入與輸出的函數關系,從而無法確定最佳參數值;如果忽略噪聲因子的影響,將會導致優化結果具有偶然性,甚至出現較大的偏差。然而,如何同時考慮模型參數的不確定性以及生產過程的噪聲因子導致的預測響應值波動問題,目前可行的研究方法還不多。

因此,針對上述的問題,本文基于統一的貝葉斯多元回歸模型框架,首先,通過蒙特卡羅模擬抽樣得到各個響應的后驗抽樣值;然后根據帕累托最優原則計算帕累托最優前沿,并獲得相應的多目標優化解集;最后,運用灰色關聯分析方法對多目標優化解集進行關聯度分析,選擇關聯度最大的作為理想的參數設計值。

1 基于貝葉斯多元回歸模型的帕累托最優前沿

1.1 貝葉斯多元回歸模型

在多響應優化設計中,假設有p個響應和q個因子效應,則多響應曲面回歸模型為

y=Bz(x)+e

(1)

式中,y為p×1的響應矩陣;B為p×q的回歸系數矩陣;z(x)為q×1的因子效應矩陣;e為p×1的隨機誤差矩陣,并且服從均值向量為0、方差-協方差矩陣為Σ的正態分布。

(2)

式中

(3)

(4)

(5)

(6)

(7)

1.2 利用貝葉斯后驗樣本的帕累托最優前沿

在數學上,一個多響應優化問題[35]可表示為

S={x|h(x)=0;g(x)≤0}

(8)

式中,h(x)為等式約束;g(x)為不等式約束;S為定義域集合。

目前,針對多響應優化問題,普遍接受的一種解釋是帕累托最優概念,即多響應優化問題不是為了求解某個唯一的最優解,而是針對響應之間的沖突關系給出一組折中解。而每一個折中解在目標空間中都對應著一個帕累托最優前沿點,這些前沿點的集合構成了帕累托最優前沿解集[30]。另外,在這個解集中,所有的點都具有帕累托最優性。

一個點x*具有帕累托最優性,即不存在這樣的點x∈S,對于所有的r∈{1,2,…,n},都有fr(x)≤fr(x*),且至少有一個為嚴格不等式[36]。

假設試驗考察的變量為x1和x2,響應為y1、y2和y3,其均為望大質量特性的響應。如果在實際研究過程中,存在望小或者望目質量特性的響應,可對數據做簡單的處理(比如加負號),將其轉化為望大質量特性的響應。根據式(7)獲得響應的貝葉斯后驗樣本抽樣值后,第k(k=1,2,…,n)個試驗點x(k)能夠達到帕累托最優的概率可以通過式(9)和式(10)近似地求得。

∈PF(i))

(9)

(10)

2 灰色關聯分析

我國學者鄧聚龍在1982年創建了灰色系統理論[37]。灰色關聯分析是灰色系統理論的重要組成部分,現在已經被廣泛應用于工程實踐中。在灰色關聯分析理論中,判斷序列之間聯系是否緊密的依據為序列曲線幾何形狀的相似程度。如果曲線幾何形狀的相似程度高,那么序列之間的聯系就大,反之就小。早期主要以鄧聚龍教授提出的灰色關聯分析模型為代表,之后許多學者在其基礎上提出多種灰色關聯分析模型。

2010年,文獻[37]在廣義灰色關聯模型的基礎上提出灰色接近關聯度模型,灰色接近關聯度用于測量序列Xi和Xj在空間中的接近關聯程度。Xi和Xj之間越接近,那么其之間的接近關聯度ρij就越大,反之就越小。以灰色接近關聯度理論為基礎,將多響應優化得到的m組帕累托解{Yi}作為比較序列{yi(1),yi(2),…,yi(n)},i=1,2,…,m,將由n個單目標各自的最優解組成的{Y0}作為參考序列{y0(1),y0(2),…,y0(n)}。令

(11)

根據式(11)可計算序列Yi和Y0的灰色接近關聯度為

(12)

式中,如果Yi和Y0長度相同,那么有式(13)成立。

(13)

在多響應優化設計中,灰色關聯分析法能夠將多響應優化問題轉化為以灰色關聯度為目標的單響應優化問題,并通過計算m組比較序列與參考序列之間的關聯度來獲得最佳因子組合和最佳響應目標值。

3 結合帕累托最優原則與灰色關聯分析的多響應優化設計

在多響應優化設計的過程中,難點之一是如何解決多個響應之間可能存在的沖突問題。目前解決這類問題的方法主要有兩種[38]:一種方法是選取其中一個響應作為目標函數,然后將剩余的響應作為約束條件,以此將多響應優化問題轉化為求解有約束的目標函數的最優解問題;另一種方法是選擇一個綜合的目標函數,然后求解該函數的最優解,由此多響應優化問題被轉化為單響應優化問題。然而,上述兩種方法各自都存在一定的缺陷。選取其中一個響應作為目標函數或者選擇一個綜合的目標函數,這兩種方法在很大程度上都受到主觀因素的影響,而且不同的目標函數往往會得到不同的優化結果。因此,利用上述兩種方法進行多響應優化設計,得到的優化結果通常是不完整的。利用帕累托最優原則的優勢,構建多響應優化算法可以彌補上述方法的不足,從而能夠客觀地得出所有符合優化目標的參數值。故本文在貝葉斯多元回歸模型統一框架下,針對在多響應優化設計中,模型參數的不確定性以及生產過程的噪聲因子導致的預測響應值波動問題,首先采用帕累托最優原則獲取一系列優化的參數設計組合,然后利用灰色關聯分析方法,進一步確定最佳的參數組合。所提方法的基本流程如圖1所示,具體步驟如下。

圖1 所提方法的多響應優化設計流程圖Fig.1 Flow chart of the proposed method for multi-response optimization design

步驟1采用中心復合設計的方法進行相關試驗,根據試驗條件,確定各個變量的試驗區域。

步驟2構建貝葉斯回歸模型統一框架,進行模型擬合優度檢驗。如果該模型通過擬合優度檢驗,則轉入步驟3,反之,需調整模型結構,轉入步驟3。

步驟4.1確定當前的抽樣次數為m,如果m≤N,其中N為抽樣總次數,則確定第m次模擬抽樣后n個試驗點的各個響應后驗抽樣值為

步驟4.2按照帕累托最優原則確定初始的帕累托最優前沿為

式中,1≤b≤i≤c≤j,此處令b=i=c=j=1。

步驟4.7得到第m次抽樣后的帕累托最優前沿PF(m)。令m=m+1,轉入步驟4.1。

圖2 計算帕累托最優前沿的算法流程圖Fig.2 Flow chart of the algorithm for calculating the Pareto optimality frontier

步驟5利用步驟4的結果結合式(9)和式(10),計算每一組參數值在整個優化過程中能夠達到帕累托最優的概率。

步驟6根據步驟5的結果,選擇概率較大的m組參數值,利用其對應的響應后驗抽樣值進行灰色關聯度分析,并按照灰色接近關聯度進行排序,以關聯度最大的為最優設計方案。在圖2的基礎上,結合前面5步的結果,灰色關聯分析的算法流程如圖3所示。

圖3 灰色關聯分析的算法流程圖Fig.3 Flow chart of grey relational analysis

步驟6.6得到ρ(s),其中s=1,2,…,m。

步驟6.7選擇ρ(s)最大的為最優設計方案,優化過程結束。

4 實例分析

4.1 實例背景

本試驗旨在通過加拿大卡爾加里大學機械與制造工程系飛秒激光微納加工中心實現某衛星芯片的微孔制造過程。該試驗方案的設計以及數據收集工作系本文第一作者在加拿大卡爾加里大學從事博士后研究工作期間,與加拿大卡爾加里大學飛秒激光微納加工中心的研究人員共同完成。該飛秒激光微納加工中心主要由飛秒激光束發生器(femtosecond laser beam generator),微加工工作站(micro-machining workstation),工控制計算機(machining control computer),測量儀器(measuring instrument)和花崗巖平臺(granite platform)組成,具體實物如圖4所示。在該試驗中,選擇直徑y1和圓度y2作為關鍵的質量特性以反映某芯片微孔的幾何特性和加工精度,其中直徑為望目質量特性,而圓度為望大質量特性。此外,圓度與微孔的面積和主軸有關,其具體的公式定義為

(14)

式中,通過圖像軟件Imagej分析微孔圖片,可以獲得式(14)中的微孔主軸和面積。影響上述兩個響應的可控因素為平均功率x1(mw)、脈沖頻率x2(Hz)和切削速度x3(mm/s),3個變量的水平和編碼設置如表1所示。

圖4 飛秒激光微納加工中心Fig.4 Femtosecond laser micro/nano-machining center

參數水平-1.682-1011.682x115.9150100150184.09x2398500650800902x30.030.060.110.160.19

本次試驗的目的是希望確定3個變量的最佳設置,從而使直徑y1的目標值維持在21的水平上,并盡可能最大化圓度y2。本文采用中心復合設計開展相關試驗,試驗數據由實驗人員通過實際的激光微鉆實驗收集,結果如表2所示。

表2 中心復合設計試驗結果

4.2 確定貝葉斯多元回歸模型統一框架

在上述試驗中,3個變量的可接受范圍均為-1.68~1.68,取值間隔為0.38,整個試驗區域包含729個試驗點。在整個優化過程中,假設回歸模型式(1)中的因子效應向量為

利用最小二乘法進行模型參數估計,計算結果如表3所示。

表3 模型擬合優度檢驗

表4 模型y1和y2兩次擬合優度對比

4.3 根據帕累托最優原則確定帕累托最優前沿

假設兩個響應的目標值分別為21和1,試驗的模擬抽樣次數N為2 000次,根據式(3)~式(7)在試驗區域內的每個點進行蒙特卡羅隨機抽樣,得到各自響應的后驗抽樣值。由于在實際操作過程中,圓度y2不存在大于1的值。因此如果該響應存在大量超過1的后驗抽樣值,需要將對應的試驗點刪除,最終試驗區域中的待測試驗點為603個。然后利用剩余試驗點的后驗抽樣結果在試驗區域內根據圖2進行多響應優化。在R優化代碼中,以矩陣的形式呈現隨機變量W和U的模擬抽樣結果。最后根據每次優化得到的帕累托最優前沿,根據式(9)和式(10)計算每個試驗點能夠達到帕累托最優的概率,其研究結果如圖5所示。

圖5 按照帕累托最優原則確定帕累托最優前沿Fig.5 Determining the Pareto optimal front based on the Pareto optimal principle

圖5中的圓點為整個試驗過程中所有分布在帕累托前沿上的試驗點,即這些點在整個優化過程中至少有一次出現在帕累托前沿上。每個圓點的大小和顏色與該點在整個優化過程中能達到帕累托最優的概率有關,點越大,顏色越深,在該點處能達到帕累托最優的可能性越大。另外,圖中的圓點大小不同,顏色深淺不一,這表明預測響應值的波動問題對優化結果產生了較大影響,這種波動問題主要來源于模型參數的不確定性以及生產過程的噪聲因子。由圖中點的分布位置可知,最有可能達到帕累托最優的試驗點主要分布在圖形的下方區域。

本例提取了整個試驗過程中前8個最有可能達到帕累托最優的試驗點,這8個試驗點對應的響應后驗抽樣值變化見圖6。

圖6 8個試驗點關于兩個響應的箱形圖Fig.6 Boxplots of two responses on the eight test points

比較分析圖6中各個試驗點針對不同響應的優化效果可知,在上述多響應優化設計的實例中,如果僅考慮響應y1上的優化效果,則第294號、329號、372號和603號試驗點給出的結果明顯優于其他試驗點。如果僅考慮響應y2上的優化效果,則以上8個點與目標值都存在一定的差距,第293號、327號、328號和372號試驗點的優化效果略高于其他試驗點。因此,根據圖6可知,除了第373號試驗點,這些可重復性較高的試驗點至少在一個響應上能夠得到較好的優化效果,這與帕累托最優的定義一致。此外,表5給出了8個試驗點各自能達到帕累托最優的概率。

分析表5的優化結果可知,在本文的多響應優化設計實例中,第372號試驗點能達到帕累托最優的可能性相對較高,而第373號試驗點的可重復性概率最低。圖6中第372號試驗點在兩個響應上都能得到較好的優化效果,而第373號試驗點關于兩個響應的優化效果相對較差,可見表5和圖6給出的結果一致,進一步說明了該優化結果是可靠的。

表5 圖6中的試驗點能夠達到帕累托最優的概率

4.4 基于灰色關聯分析的多響應優化

在上述試驗的優化過程中,運用本文所提方法獲得了8組優化結果。此處運用灰色關聯分析方法,對表5中的8組結果進行排序。為此,將本例中各個響應的目標值{21.000,1.000}作為參考序列,根據圖3計算每組解對應的響應后驗抽樣值與參考序列之間的灰色接近關聯度,計算其平均值作為該組解與參考序列的灰色接近關聯度,按照此關聯度對8組解進行排序,以關聯度最大的作為最優設計方案。灰色關聯分析法的研究結果如表6所示。

表6 灰色關聯分析結果

結合表6的分析結果可知,第372號、329號和第294號試驗點的灰色關聯度非常接近,這3個也是可重復性概率相對較高的試驗點。其中,第372號試驗點的關聯度略大于其余兩個試驗點,應為本次試驗的最優結果,其對應的兩個響應值分別為20.886 9和0.968 1。此外第603號和第373號試驗點的關聯度較接近,而第327號、第328號和第293號試驗點的關聯度差距不大。雖然前者在響應y2上的優化效果略低于后者,但是在響應y1上的優化效果明顯優于后者。所以在表5中排名靠后的第603號和第373號試驗點的關聯度略大于第327號、第328號和第293號試驗點。可見,針對預測響應值存在波動的多響應優化問題,運用本文所提的方法,能夠從可重復性和關聯度兩個方面獲得較為滿意的結果。

5 結 論

在多響應優化設計中,往往需要考慮預測響應值波動對優化結果的影響。本文在貝葉斯回歸模型的統一框架下,結合帕累托最優策略和灰色關聯分析方法,提出了一種多響應優化設計的新方法。該方法在多響應優化設計過程中同時考慮了模型參數不確定性和生產過程中噪聲因子的影響,并結合貝葉斯后驗概率方法計算試驗點能夠達到帕累托最優的概率,然后利用灰色關聯分析方法進一步確定理想的參數值。另外,通過實例證明該方法對解決預測響應值存在波動的多響應優化問題是有效的。

需要特別指出的是,本文中試驗區域的選擇具有主觀性,該案例僅有3個變量,2個響應,試驗區域較小,利用蒙特卡羅模擬抽樣方法和帕累托最優策略能夠快速地得到帕累托最優前沿。如何在變量和響應個數增加時選擇合適的試驗區域,有待進一步地深入研究。此外,本文的研究基于貝葉斯多元回歸模型的統一框架,在模擬優化過程中,假設各個響應的模型結構完全一致,如何利用SUR模型進行多響應優化設計,以考慮各個響應模型的結構對多響應優化結果的影響,這是今后需要研究的課題之一。另外,盡管SUR模型能夠更精確地擬合響應模型,但是在SUR模型的貝葉斯建模過程中其模擬響應的抽樣速度非常緩慢,這也是未來研究過程中有待解決的關鍵科學問題之一。

猜你喜歡
優化方法模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
主站蜘蛛池模板: 国产亚洲欧美在线专区| 丁香婷婷综合激情| 国产第一页免费浮力影院| 日韩精品免费一线在线观看| 白丝美女办公室高潮喷水视频 | 色国产视频| 91精品国产综合久久香蕉922| 国产杨幂丝袜av在线播放| 精品国产成人av免费| 国产人成乱码视频免费观看| 亚洲成AV人手机在线观看网站| 欧美黄色网站在线看| 日韩经典精品无码一区二区| 国产欧美日韩在线一区| 五月天久久婷婷| 国产激情在线视频| 亚洲AV无码久久天堂| 免费看美女毛片| 午夜精品久久久久久久无码软件| 天天躁夜夜躁狠狠躁躁88| 韩国v欧美v亚洲v日本v| 18禁黄无遮挡网站| 四虎亚洲国产成人久久精品| 国产一区亚洲一区| h视频在线播放| 亚洲第一黄片大全| 国产精品久久精品| 中文字幕2区| 久久精品人妻中文视频| 国产网站黄| 成人免费黄色小视频| 亚洲区欧美区| 国产精品浪潮Av| 最新国产精品鲁鲁免费视频| 亚洲色图欧美在线| 亚洲第一视频网站| a毛片在线免费观看| 最新亚洲av女人的天堂| 青青操视频在线| 欧美日韩另类国产| 久久99热这里只有精品免费看| 97国产在线视频| 成人av手机在线观看| 一级毛片无毒不卡直接观看| 国产免费观看av大片的网站| 国产自在自线午夜精品视频| 国产又粗又猛又爽| 久久香蕉国产线看观看亚洲片| 亚洲永久精品ww47国产| 欧美午夜在线播放| 国产成人AV综合久久| 国产三区二区| 亚洲一区精品视频在线| 激情综合网激情综合| 亚洲第一黄色网| 四虎亚洲国产成人久久精品| 亚洲品质国产精品无码| 国产男女免费完整版视频| 不卡的在线视频免费观看| 亚洲一级毛片免费看| 波多野结衣久久高清免费| 九九九精品成人免费视频7| jizz国产视频| 色婷婷成人| 国产精品入口麻豆| 亚洲天堂视频网站| 亚洲无码37.| 国产精品999在线| 国产福利小视频在线播放观看| 免费高清毛片| 丁香亚洲综合五月天婷婷| 99国产精品国产| 又粗又硬又大又爽免费视频播放| 国产一级小视频| 亚洲日韩国产精品无码专区| 欧美一区国产| 国产成人乱无码视频| 18禁影院亚洲专区| 中文精品久久久久国产网址| 97国产在线播放| 亚洲男人天堂网址| 美女无遮挡免费网站|