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

一種新型Powell粒子群算法同步綜合換熱網絡

2016-10-22 07:14:31張春偉崔國民
化工進展 2016年10期
關鍵詞:優化

張春偉,崔國民

(上海理工大學新能源科學與工程研究所,上海 200093)

一種新型Powell粒子群算法同步綜合換熱網絡

張春偉,崔國民

(上海理工大學新能源科學與工程研究所,上海 200093)

針對換熱網絡同步綜合方法的不足,本文提出了一種新型Powell粒子群算法,具有傳統確定性方法的高精度以及啟發式方法的高效率。同時針對群體智能算法優化換熱網絡問題時存在的不足,提出了云記憶體和個體對立策略,有效地避免算法發生早熟現象,擴大搜索范圍。為處理整型變量而提出的兩條整型變量優化策略與 Powell粒子群算法結合,實現了連續變量與整型變量的同步優化。最后,選取兩個經典算例驗證算法的性能,均獲得了優于文獻的結果,表明算法能夠找到更優的換熱網絡結構,是一種處理混合整數非線性問題的有效方法。

換熱網絡綜合;整型變量優化;Powell法;群體智能

換熱網絡綜合(heat exchanger network synthesis,HENS)是過程系統工程的一個重要領域,其目的是發現一個能量回收最大或者年綜合費用(包括投資費用和運行費用)最小的換熱網絡結構,進而提高能量利用率和經濟性。換熱網絡綜合問題的復雜性很大程度源于其換熱器的組合本質,FURMAN等[1]證明其為NP-難問題,因此即使是小規模的換熱網絡問題也很難證實得到全局最優解。

換熱網絡綜合方法可分為分步綜合方法和同步綜合方法兩類[2]。夾點技術法[3]由于其操作簡便、觀點明確,是一種在工程上應用較多的分步綜合方法。但由于不能很好地權衡投資費用與能量回收,只能獲得次優解。同步綜合方法[4]的數學模型屬于混合整數非線性規劃范疇(mixed integer non-linear programming,MINLP),而基于同步綜合模型的全局優化方法可進一步分為確定性方法和啟發式方法[5]。確定性方法如分支定界法[6]、外部逼近法[7]等,但其在處理大規模換熱網絡問題時效率低下。啟發式方法已在換熱網絡問題中得到了一定的應用,如遺傳算法[8]、微分進化算法[9]、粒子群算法[10]等。此類算法通過仿生自然過程具有較強的全局搜索能力,但易受種群多樣性等因素的影響,發生早熟等現象,進而陷入局部最優解。此外,由于不需要目標函數的梯度等信息,所以其求解精度有限。

在優化換熱網絡這類非凸、非線性嚴重的MINLP問題時,效率和精度一直認為是評價算法性能的兩個重要指標,而常用的同步綜合方法受限于其優化機理,很難兼顧兩者。傳統的確定性方法對初始點的依賴性較大,只能得到局部區域內的最小解。但依賴于梯度或方向的優化特性,使其具有較高的精度。所以在一些研究[11-13]中,此類算法會作為啟發式方法的局部搜索策略,即在前者的搜索機理不變的情況下,增加算法的精度。但此種結合方式未充分發揮兩種算法的優點,在優化換熱網絡時,會使算法的計算時間無限增長并忽略大量的較優解。

鑒于此,本文提出了一種新型 Powell粒子群算法(Powell particle swarm optimization algorithm,PPSO),將傳統確定性方法的高精度、收斂速度快的特性與啟發式方法的全局搜索能力相結合,同步綜合換熱網絡。在算法中,每個粒子均具有Powell法的局部優化能力,通過構造新的啟發式準則將所有粒子組織起來,使其具備全局搜索能力以及移動特性。并針對換熱網絡問題以及算法特性,提出了云記憶體和個體對立策略,擴大算法的搜索范圍,避免早熟現象的發生。同時為處理問題中的整型變量,提出兩條相關的整型變量處理策略。最后,通過兩個經典算例對算法的性能進行了驗證。

1 換熱網絡數學模型

1.1 換熱網絡問題數學描述

換熱網絡綜合問題可表述如下:現有 NH股熱流體、NC股冷流體,分別需要冷卻、加熱到相應的目標溫度。在冷、熱流體之間設置多個換熱器,實現能量回收。當某一股流體未達到目標溫度時,為其匹配冷或熱公用工程。過程流體和公用工程的熱容流率、進出口溫度以及換熱器換熱系數均已知。以年綜合費用為優化目標,包括運行費用和投資費用兩部分,運行費用為消耗冷、熱公用工程時產生的費用,投資費用為設置換熱器時產生的費用,可分為面積費用和固定投資費用。

對于有分流的換熱網絡模型而言,雖然有時能夠得到相對無分流模型費用值更低的換熱網絡設計,但由于其增加了問題的復雜性,而且沒有計算由于流體分流引起的管路、閥門等造成的額外投資費用,所以在實際工程中的可接受度不高。鑒于此,本文采用GROSSMANN等[4]提出的無分流分級超結構模型,其中,換熱網絡的級數為NS=max(NH,NC),換熱器的最大個數為NK=NS×NH×NC。現以2股熱流體和2股冷流體為例表述分級超結構模型,如圖1所示。

圖1 換熱網絡無分流的分級超結構

1.2 優化的目標函數

針對上述換熱網絡模型,以年最小年綜合費用為優化目標,其數學函數為式(1)。

式中,B為0-1整型變量,表示換熱器有無。當換熱器存在時,B=1,反之,則B=0;CCU、CHU分別為冷、熱公用工程的費用系數;CF為設置換熱器時的固定投資費用;CE為面積費用系數;Z為面積費用指數;QHU,j、QCU,i分別為熱公用工程與冷流體之間的換熱量和冷公用工程與熱流體之間的換熱量;AHU,j、ACU,i分別為其相應換熱器的面積;Ai,j,k為冷熱流體之間換熱器面積。本文以單個換熱器的換熱量Qi,j,k為優化變量,各換熱器均采用逆流傳熱方式,如式(2)~式(4)。

1.3 約束條件

2 Powell粒子群算法

Powell粒子具有找到所在區域局部最優的能力,所以粒子更容易發生聚集。鑒于此,本算法的中心思想是使粒子最大程度地遍歷搜索空間同時避免發生早熟現象。Powell法一種求解無約束最優化問題的直接搜索法[14],具有收斂速度快、求解精度高、無需計算導數等優點,較其他傳統的局部優化方法更適合換熱網絡綜合問題。Powell法的優化換熱網絡問題時的計算步驟如下所示。

2.1 Powell法計算步驟

Step 1 確定變量維數N,即粒子結構的所包含的換熱器個數,讀取當前位置信息Q0。設置收斂精度ε1。給定一組線性無關的方向 Di(i=1,2,…,N),Di取個坐標軸的N個方向,即N階單位矩陣。其中N的最大取值為NK。

Step 2 從初始點 Q0出發依次沿方向 Di(i=1,2,…,N)進行一維搜索,確定每次迭代的步長 ,得到Q1,Q2,…,QN,見式(32)、式(33)。

Step 3 判斷迭代計算是否結束:若滿足式(34),則得到解QN,計算結束;否則轉Step 4。

Step 4 計算最速上升方向上函數 F(Q)的變化,見式(35)。

Step 5 引進第(N+1)個搜索方向和新的點Qt,見式(36)、式(37)。并計算F(Qt)。

Step 6 方向替換判斷。

①若滿足

則將 QNh作為新的初始點,沿原方向搜索,即轉Step 2。

③若滿足

則將 QN作為新的初始點,沿原方向搜索,即轉Step 2。

③若以上兩條件均不滿足,則轉Step 7。

Step 7 以QN作為起始點,沿方向DN+1進行一維搜索,并得到此方向上的極小值點QN+1。將方向Dibig用新方向DNh+1替換,產生一組新的方向Di(i=1,2,…,N),以QN+1作為新初始點,轉Step 2。

2.2 粒子的更新公式

根據算法的中心思想以及群體智能算法,本文提出了如式(40)的Powell粒子更新公式。

式中,r1和r2為介于(0,1)之間的偽隨機數;1β和2β為非負常數;Vn,max為粒子的最大飛行速度;為粒子自身存儲的最優位置;為粒子鄰居的最優位置。公式第一項為隨機擾動項,調整粒子的飛行速度和方向。當粒子的飛行速度較大時,其本身的隨機搜索能力較弱,當粒子本身速度較小時,則產生一個相對較大的隨機向量。由于算法是在保證種群進行信息交流的前提下,最大程度地遍歷搜索空間,所以此項可以避免算法發生早熟現象。第二項與第三項分別為個體向自身最優與鄰居最優的飛行的方向向量。即Powell粒子會根據隨機搜索項、自身最優以及群體最優位置調整下一次的飛行方向并記憶搜索到的最優位置。

2.3 馮諾依曼拓撲結構

拓撲結構影響著算法的收斂速度以及精度等,其總體上可分為全局拓撲和局部拓撲兩種。全局拓撲結構雖然收斂速度快,但易于陷入局部最優解。鑒于此,本算法采用局部拓撲中的馮諾依曼拓撲結構,如圖2所示。每個Powell粒子只有4個鄰居,相應的為鄰居中的最佳位置。當一個粒子找到較好解時,只影響其周圍的4個個體,這樣可以較好地維持種群多樣性。使得種群中的粒子向不同的最優位置靠攏,不易陷入局部最優。計算結果表明,此種拓撲結構對于復雜的問題具有良好的收斂性能,適用于解決換熱網絡綜合問題。

圖2 馮諾依曼拓撲結構

2.4 云記憶體

在群體智能算法中,個體的記憶能力是一個很重要的部分,能夠指導個體搜索到更優的位置。如果算法采用常規的信息存儲方式,即個體只記憶自身搜索到的最佳位置,可能會導致種群搜索到的一些較優解被忽略。例如,記錄著種群當前最優位置的Powell粒子的記憶很難更新,除非發現了新的種群最優解,那么差于當前最優位置卻優于其他個體記憶位置的解則會被算法遺失。

鑒于此,本文將相互獨立的個體記憶功能統一組合為整個種群的記憶功能,即云記憶體,其內部保存的信息個數等于種群中的個體數,并且所有信息已按照適應度值的大小進行排序。當一個新解產生后,與種群記憶體中適應度值最差的位置信息進行比較,若優于此解,則將新解與記憶體中的所有信息再次進行排序,進而剔除最差的位置信息;若差于此解,則記憶體保持不變,算法繼續進行搜索。根據云記憶體的思想,個體與其存儲的最優位置信息不再有直接的聯系,即其存儲的信息可能為自身搜索到的,也可能為其他個體搜索到后存儲在當前個體上。以圖3所示為例說明云記憶體存儲方式,假設種群中共有4個Powell粒子,那么云記憶體中則包含4個最優位置信息,并分別存儲在不同的粒子上。通過云記憶體,算法能夠克服常規信息存儲方式的不足。此外,本文的粒子更新公式基于自身與鄰居最優位置,所以云記憶體也會增加算法的隨機性,擴大搜索范圍。

圖3 種群云記憶體

2.5 個體對立策略

Powell粒子能夠得到局部區域的極小值,所以當某些粒子處在一個相同的區域時,可能會發生聚集,弱化算法的性能。所以為了提高算法的搜索效率同時擴大搜索范圍,提出個體對立策略,處理陷入局部最優的粒子。粒子狀態的判定公式如式(41)、式(42)所示。

Smove表示Powell粒子的當前移動距離,當其小于一個設定的值時,認為此粒子陷入局部最優區域,此時通過個體對立公式變換粒子的位置,如式(43)。

式中,Qn,K為Powell粒子的當前位置;為對立后產生的新位置;為粒子n的搜索區間上限,即為換熱器所在流體的最大換熱潛能。由于算法能夠不斷地跳出局部最優解,所以具有很強的“爬山”能力。

3 整型變量優化策略

隨著換熱網絡問題規模的日益增大,導致整型變量急劇增加,可行結構數量呈幾何級增長,所以同步綜合換熱網絡問題時,有效的整型變量優化技術不可或缺。鑒于此,本文提出了兩條相關的整型變量優化策略。為準確表示換熱網絡結構,此處采用一種數字序列表達方式,兩者的轉換公如式(44)所示。由其可知,當換熱器不存在時,序列中對應的數字為零,所以為表述簡潔,序列中只包含存在的換熱器,零則不予表述。詳細轉換過程如圖 4所示。

圖4 4股流換熱網絡結構及其序列表示

3.1 整型變量判斷策略

在算法初始化階段,設定個體的初始解均為全結構,即所有換熱器都存在。通過測試發現,計算過程中,粒子某一維的換熱量Q近似為零。從投資費用方面分析,當換熱器上的熱負荷小于一個特定的值時,可以認為此換熱器對于減少年綜合費用值是不利的,即這一維表示的換熱器應該消去。但由于精度的影響,優化變量Q不可能完全為零。鑒于此,本文提出了一條整型變量判斷策略,其形式如式(45)所示。

式中,Qmin為設定的換熱器最小換熱量,當Powell粒子某一維的變量小于此值時,令序列中的相應位置的整數為零,即消去對應的換熱器。

3.2 整型變量進化策略

整型變量進化策略是對整型變量判斷策略的進一步延伸和拓展。但從費用函數角度分析,過多的換熱器是不利于費用值降低的,所以本文以整型變量判斷策略為基礎,提出了整型變量進化策略。即Powell粒子進行位置更新時,對其每一維所表示的換熱器進行隨機消去操作。詳細操作如下:設定消去概率 pe,每一維變量進行更新時,產生隨機數rand,當rand≤pe時,消去序列中整數表示的換熱器。對于新產生的結構,本文通過構建一個接受概率函數pr,其形式如式(46)所示。

式中, Fold和Fnew分別為原結構和新結構對應的目標函數值,α=0.95為縮減因子,p0=0.2為pr的初始值。此處應注意的是,當換熱器被消去后,若新結構的費用值降低,則此新解作為最佳的概率相應增大,種群中其他個體傾向于消去此換熱器;若其費用值升高,則此新解作為最佳的概率相應減小,所以其他個體消去此換熱器的概率也相應減少,而此個體在下次迭代過程中又可以通過學習其他個體,重生已消去的換熱器。所以整型變量進化策略是一種以費用值為導向的進化策略,整型變量判斷策略與其類似。

3.3 算法流程

Powell粒子群算法與整型變量優化策略相結合,能夠實現換熱網絡的連續變量與整型變量的同步優化。其主要步驟為:隨機生成 Np個 Powell粒子,包括初始位置向量Qn,0和速度向量V,其中,每個粒子均代表了一組潛在解。粒子在整個搜索空間內根據隨機搜索特性、自身最優和鄰居最優位置進行位置更新,同時執行整型變量優化策略,如圖5所示。

圖5 算法流程圖

4 算例驗證與分析

為驗證Powell粒子群算法的性能,本文采取兩個典型算例對其進行驗證。計算環境為Win7系統下Fortran(Compaq Visual Fortran 6),計算機參數為Intel(R) Xeon(R) CPU E5-2670 v2 2.5GHz 32GB RAM。

4.1 算例一

算例一取自文獻[15],包含6股熱流體與5股冷流體,算例的相關計算參數如表 1所示。CASTILLO 等[15]所得結構對應的費用值為141555$/a。SILVA 等[10]采用一種內、外層均為粒子群的雙層算法同步綜合換熱網絡,其獲得的有分流分級超結構對應的年綜合費用139777$/a。但在本文的計算中發現,在H4C3流股上的換熱器的溫度出現交叉,即熱流體 H4的出口溫度比冷流體 C3的入口溫度低,如文獻[10]中的圖 6所示,因此認為他們的結構與費用值不符合,所以其計算結果不在本文的比較范圍內。

表1 算例一參數

采用Powell粒子群算法優化算例一,得到圖6所示的結構,變量單位為 kW,其費用值為140084$/a。通過分析表1中的數據可以發現,其熱容流率數值差異很大,所以導致此算例的非線性很強,過程流體參數的細微改變,都會導致費用值的急劇變化。所得結果與文獻的對比情況如表2所示,由其可知,兩者的換熱單元數與消耗的公用工程量近似相同,但根據參考文獻[15]的結構發現,兩者的換熱網絡結構并不同,即本文獲得了優于文獻的結果,證明了算法的高精度。

圖6 算法優化結果(年綜合費用:140084$/a)

表2 算例一結果比較

4.2 算例二

算例二取自文獻[16],是由 10股熱流體與 10股冷流體組成大規模的換熱網絡問題,算例的相關參數如表3所示。XIAO等[16]采用一種傳統的溫焓圖法設計多股流換熱網絡,其所得結果為1827772$/a。LUO 等[17]首次將此算例應用于兩股流換熱網絡設計,并采用一種由遺傳算法、模擬退火算法等組成的混合算法優化算例二,其獲得的有分流分級超結構的費用值為1753271$/a。LAUKKANEN等[18]提出了一種對過程流體進行分組的雙層算法來減少換熱網絡綜合問題的復雜性,其獲得的結構費用值為1811.9k$/a。

表3 算例二參數

采用Powell粒子群算法優化算例二,得到圖7所示的結構,變量單位為 kW,其費用值為1745145$/a,與文獻的結果對比如表4所示。可以發現,本文獲得了相對于有分流分級超結構更優的結果,證明了算法處理大規模換熱網絡綜合問題的高效性。計算時的相關參數如表5所示,所有的參數均為多次計算所取的經驗值。

圖7 算法優化結果(年綜合費用:1745145$/a)

表4 算例二結果比較

表5 計算相關參數

5 結 論

本文提出了一種新型Powell粒子群算法,粒子的Powell特性使算法具有較高的搜索精度,而針對換熱網絡綜合問題構建的啟發式準則使算法具有較強的全局搜索能力,能夠同時兼顧精度和效率。為優化換熱網絡結構而提出兩條整型變量優化策略能夠與Powell粒子群算法很好的契合,可以有效地處理問題中的整型變量。采用兩個經典算例對算法進行驗證,均獲得了相對文獻更優的結果,證明了算法的性能。

符 號 說 明

A——換熱器面積,m2

B——取值為0或1的邏輯變量

C——費用計算系數

NC——冷流股數

NH——熱流股數

NS——換熱網絡級數

NK——最大換熱器數目

Np——種群個數

Q——換熱器換熱量,kW

V——粒子的速度

rand——取值介于(0,1)的隨機數

T——溫度,℃

Z——面積費用指數

上角標

in——換熱器入口

out——換熱器出口

下角標

c —— 冷流體

h —— 熱流體

i —— 熱流股標號

j —— 冷流股編號

k —— 級數編號

CU —— 冷公用工程

HU —— 熱公用工程

[1] FURMAN K C,SAHINIDIS N V. Computational complexity of heat exchanger network synthesis[J]. Computers & Chemical Engineering,2001,25(9):1371-1390.

[2] FURMAN K C,SAHINIDIS N V. A critical review and annotated bibliography for heat exchanger network synthesis in the 20th century[J]. Industrial & Engineering Chemistry Research,2002,41(10):2335-2370.

[3] LINNHOFF B,HINDMARSH E. The pinch design method for heat exchanger networks[J]. Chemical Engineering Science,1983,38(5):745-763.

[4] YEE T F,GROSSMANN I E. Simultaneous optimization models for heat integration——Ⅱ. Heat exchanger network synthesis[J]. Computers & Chemical Engineering,1990,14(10):1165-1184.

[5] CHOI S H,MANOUSIOUTHAKIS V. Global optimization methods for chemical process design:deterministic and stochastic approaches[J]. Korean Journal of Chemical Engineering,2002,19(2):227-232.

[6] GROSSMANN I E,SARGENT R W H. Optimum design of multipurpose chemical plants[J]. Industrial & Engineering Chemistry Process Design and Development,1979,18(2):343-348.

[7] ZAMORA J M,GROSSMANN I E. A global MINLP optimization algorithm for the synthesis of heat exchanger networks with no stream splits[J]. Computers & Chemical Engineering,1998,22(3):367-384.

[8] DIPAMA J,TEYSSEDOU A,SORIN M. Synthesis of heat exchanger networks using genetic algorithms [J]. Applied Thermal Engineering, 2008,28(14):1763-1773.

[9] YERRAMSETTY K M,MURTY C V S. Synthesis of cost-optimal heat exchanger networks using differential evolution [J]. Computers & Chemical Engineering,2008,32(8):1861-1876.

[10] SILVA A P,RAVAGNANI M A S S,BISCAIA Jr E C,et al. Optimal heat exchanger network synthesis using particle swarm optimization[J]. Optimization and Engineering,2010,11(3):459-470.

[11] 張安玲,王中. 一種混合粒子群優化算法的研究[J]. 計算機工程與應用,2011,47(31):27-29.

[12] 吳建輝,章兢,陳紅安. 融合Powell搜索法的粒子群優化算法[J].控制與決策,2012,27(3):343-348,354.

[13] 羅平,姚立海,楊仕友,等. 基于移動最小二乘法和粒子群算法的優化算法[J]. 浙江大學學報:工學版,2006,40(9):1482-1485.

[14] POWELL M J D. An efficient method for finding the minimum of a function of several variables without calculating derivatives [J]. The Computer Journal,1964,7(2):155-162.

[15] CASTILLO E,ACEVEDO L,REVERBERI A P. Cleaner production of nitric acid by heat transfer optimization: a case study [J]. Chemical and Biochemical Engineering Quarterly,1998,12(3):157-165.

[16] XIAO Wu,DONG Hongguang,LI Xinqing,et al. Synthesis of large-scale multistream heat exchanger networks based on stream pseudo temperature [J]. Chinese Journal of Chemical Engineering, 2006,14(5):574-583.

[17] LUO X,WEN Q Y,FIEG G. A hybrid genetic algorithm for synthesis of heat exchanger networks [J]. Computers & Chemical Engineering, 2009,33(6):1169-1181.

[18] LAUKKANEN T,FOGELHOLM C J. A bi-level optimization method for simultaneous synthesis of medium-scale heat exchanger networks based on grouping of process streams[J]. Computers & Chemical Engineering,2011,35(11):2389-2400.

A novel Powell particle swarm optimization algorithm for simultaneous synthesis of heat exchanger networks

ZHANG Chunwei,CUI Guomin
(Research Institute of New Energy Science and Technology,University of Shanghai for Science and Technology,Shanghai 200093,China)

Due to the defect of simultaneous methods for heat exchanger networks synthesis,a novel Powell particle swarm optimization(PPSO) algorithm was proposed,which has both high precision of the deterministic methods and high efficiency of the stochastic algorithms. For overcoming disadvantages of stochastic algorithms,the cloud memory and the opposite strategy of individualswere proposed,which can avoid premature convergence and expand the search space. In addition,two optimizing strategies of integer variables were combined with PPSO algorithm in order to simultaneously optimize continuous and integer variables. The presented approach was tested on two typical benchmark problems. The obtained solutions are better than that published in the literature. Results showed that the presented algorithm can find better designs,which is conductive to cost saving in industrial production.

heat exchanger networks synthesis;integer variables;Powell method;swarm intelligence

TK 124

A

1000-6613(2016)10-3092-09

10.16085/j.issn.1000-6613.2016.10.012

2016-01-27;修改稿日期:2016-02-23。

國家自然科學基金項目(51176125)。

張春偉(1992—),男,碩士研究生,從事過程系統優化研究。聯系人:崔國民,教授,博士生導師。E-mail cgm1226@163.com。

猜你喜歡
優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
PEMFC流道的多目標優化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
圍繞“地、業、人”優化產業扶貧
今日農業(2020年16期)2020-12-14 15:04:59
事業單位中固定資產會計處理的優化
消費導刊(2018年8期)2018-05-25 13:20:08
4K HDR性能大幅度優化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 免费A∨中文乱码专区| 色悠久久久久久久综合网伊人| 欧美精品成人| 国产精品国产三级国产专业不| 黄色网站不卡无码| 国内丰满少妇猛烈精品播| 亚洲成年人网| 国产男女XX00免费观看| 日韩精品亚洲人旧成在线| 亚洲无码在线午夜电影| 午夜啪啪福利| 在线亚洲天堂| 欧美一级在线看| 少妇精品网站| 69综合网| 精品国产www| 欧美午夜视频| 色网站在线免费观看| h视频在线播放| 欧美日韩亚洲综合在线观看| 怡春院欧美一区二区三区免费| 亚洲天堂视频网站| 天天综合网亚洲网站| 亚洲综合经典在线一区二区| 国产在线八区| 99re精彩视频| 日韩精品一区二区三区大桥未久| 美女被操黄色视频网站| 亚洲欧美自拍中文| 国产熟女一级毛片| 日本高清在线看免费观看| 精品一区二区三区自慰喷水| 国产成人AV男人的天堂| 99视频在线免费| 国产黄在线免费观看| 国产成人综合在线视频| 中文字幕永久视频| aⅴ免费在线观看| 毛片在线播放a| 欧美精品亚洲日韩a| 日本精品一在线观看视频| 丰满人妻久久中文字幕| 亚洲中文字幕无码爆乳| 直接黄91麻豆网站| 国产高清国内精品福利| 久久77777| 国产高清无码麻豆精品| 麻豆精品在线| 国产精品浪潮Av| 免费看a级毛片| 欧美人人干| 黄色成年视频| 国产在线视频导航| 成人在线欧美| 久草视频精品| 亚洲人成亚洲精品| 大乳丰满人妻中文字幕日本| 岛国精品一区免费视频在线观看| 免费无码在线观看| 欧美色视频日本| 久久精品欧美一区二区| 亚洲香蕉在线| 精品少妇人妻无码久久| 免费无码AV片在线观看国产| 午夜国产精品视频| 精品视频一区在线观看| 日韩天堂网| 国产剧情无码视频在线观看| 波多野结衣亚洲一区| 亚洲福利一区二区三区| 香蕉久人久人青草青草| 亚洲福利一区二区三区| 亚洲国产精品一区二区高清无码久久| 91精品国产综合久久不国产大片| 萌白酱国产一区二区| 国产美女免费网站| 久久国产香蕉| 国产特级毛片aaaaaa| 国产精品欧美亚洲韩国日本不卡| 九一九色国产| 中文字幕亚洲精品2页| 高清无码不卡视频|