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

采用中心平衡優化的表冷系統預測控制研究

2023-04-11 12:36:38盧志敏嚴德龍王華秋
關鍵詞:閥門優化模型

盧志敏,饒 偉,江 琳,嚴德龍,王華秋,向 力

(1.龍巖煙草工業有限責任公司, 福建 龍巖 364021;2.重慶理工大學 兩江人工智能學院, 重慶 401135;3.重慶太和空調自控有限公司, 重慶 400030)

表冷器溫濕度協同優化一向是暖通空調領域研究的熱門問題[1]。表冷器溫濕度協同控制是一個十分復雜的過程,表現在強耦合和干擾因素多,其過程控制是非線性、大時滯的多輸入多輸出系統[2],當使用PID控制、自適應控制、模糊控制等一些經典的控制方法時,若過程中存在某些外界或者內部的干擾時,系統的控制效果會比較難達到理想的目標[3]。預測控制由3個部分組成,先對控制對象建立預測模型,觀測其變化規律,然后利用滾動優化確定控制量,最后利用預報誤差校正預測模型[4-5]。盡管大多數預測模型可以逼近復雜的非線性函數[6],但對于空調表冷器系統,仍存在精度不高的問題[7]。

極限學習機(extreme learning machine,ELM)結構簡單、計算速度快、擬合能力強,是一種非常適合的非線性預測模型,已經被許多學者應用到工業控制過程的軟測量和軟計算中[8]。由于ELM預測模型是一種黑箱預測,因此滾動優化過程就不能采用傳統的求導方法執行,所以本文通過設置合理的目標適應度,采用帶約束的新型群智能優化算法求解最優的控制量。平衡優化算法(equilibriumoptimizer,EO)同樣具有計算簡單、收斂速度快、全局優化等優點,適合于工業過程的控制量尋優[9-10],因此本文將EO運用于表冷器預測控制系統的滾動優化環節。

針對極限學習機隱含層節點數不適當會影響預測精度和平衡優化算法容易陷入局部極小的問題,本文用斐波拉契搜索最優的隱含層節點數,采用K中心聚類算法產生候選解集,從而提出斐波拉契極限學習機(FELM)和中心平衡優化(KEO)算法,接著將FELM用于表冷系統的預測模型中,預測表冷系統下一時刻的輸出,并將KEO用于該系統的滾動優化環節中,以確定當前時刻的控制量,結合原有的誤差反饋機制,校正FELM預測模型,這樣就構成了完整的空調表冷系統的預測控制模型。

1 空調表冷系統

空調表冷系統用于調節卷煙廠工藝車間的室內溫度和室內相對濕度,通常的做法是將新回風混合,然后冷卻到較低的露點溫度,然后進行二次加熱或加濕調節室內溫度和室內相對濕度達到設定值。雖然這種做法便于操作,但是需要冷卻后再加熱加濕,存在二次加熱和加濕的問題,從而形成冷熱交替現象導致能耗巨大。

將現有表冷器進行改造,實現表冷分層以達到除濕降溫的協同控制。改裝后的分層表冷器,下層6組作為主表冷器用來除濕,上層4組作為副表冷器用來降溫。先調節主表冷電動閥,調整除濕量,當除濕滿足要求后,調節電動三通閥閥門開度來控制回水比例,主表冷器利用在除濕過程中的部分冷媒水調節溫度,然后調節副表冷電動閥閥門開度來補充冷水量,得到主副表冷器的混合冷水,達到與使用2種冷媒源相同的效果,從而實現溫濕度協同控制。

通過在表冷器上使用溫濕協同控制裝置,如圖1所示。主表冷電動閥通過控制冷凍水流量控制主表冷器除濕量。溫濕協同裝置在主表冷電動閥后設置電動可控三通閥,利用主表冷器回水2來補充副表冷器的冷凍水,兩者進行混合后一起進入副表冷器,增加副表冷器的冷水流量,同時又有效利用了主表冷器的部分回水,達到了節能的作用。副表冷電動閥則控制著進入副表冷器冷凍水,由于有主表冷器回水2的補充,副表冷的供水量也會有所下降,通過改變冷水溫度改變副表冷器的換熱量。這樣就實現了溫濕度協同控制,消除再熱現象,達到節能的目的。

圖1 表冷器溫濕協同控制裝置

在不同的季節交替運轉工況下,使用FELM建立空調表冷系統預測模型,使用改進的KEO算法滾動優化得到最優的閥門開度控制量,提出一種模式預測控制方法對分層表冷器協同控制過程中的溫濕度進行控制。

2 斐波那契ELM預測模型

2.1 ELM預測模型

極限學習機(extreme learning machine,ELM)的結構類似于三層BP網絡,但是計算方式完全不同,輸入層到隱含層之間的連接權值可以隨機確定,而不需要用反向傳遞的誤差來調節了,隱含層到輸出層之間連接權值通過求解一個包含輸出目標的矩陣獲得。由于沒有了反向傳播修改權值,而且考慮了輸出目標,使得ELM能夠快速準確地逼近任意非線性函數[11]。考慮到預測模型關系到控制器的穩定性和準確性[12],建立了基于極限學習機的表冷系統預測模型。

設計的多輸入多輸出的有約束多變量非線性系統可表示為:

y1(k+d)=f1(o1(k),o2(k),u1(k),u2(k),u3(k))

y2(k+d)=f2(o1(k),o2(k),u1(k),u2(k),u3(k))

(1)

約束條件:

20≤o1(k)≤30;50≤o2(k)≤70;

0

其中:o(k)=(o1(k),o2(k))為第k時刻環境溫濕度,作為一個外部條件輸入到預測模型中;u(k)=(u1(k),u2(k),u3(k))為第k時刻輸入向量,也就是需要優化的控制量,代表3個閥門的開度,范圍是0~1;y(k)=(y1(k),y2(k))為第k個時刻的輸出向量,也就是室內目標溫濕度;d為閥門調節后室內溫濕度的滯后時間。式(1)中的約束條件限定了控制量和外部條件的范圍。

ELM預測模型的多輸入多輸出結構如圖2所示。基于ELM的表冷系統預測模型方法的主要思想就是將主表冷閥開度u1、副表冷閥開度u2和電動三通閥開度u3,以及室外溫度o1和室外相對濕度o2作為ELM的輸入向量,將室內平均濕度y1和室內平均溫度y2作為ELM的輸出值,通過時間序列數據進行訓練,建立了ELM的5輸入2輸出表冷系統預測模型。

圖2 ELM預測模型結構

2.2 極限學習機的計算過程

極限學習機是一種特殊的前饋神經網絡[13],訓練過程中,隱含層的權值和偏置往往是隨機產生或者人為給定,不需要更新,通過計算求出輸出層的權值,即完成訓練過程[14]。步驟如下:

1) 獲取數據。數據預處理,一般是采用極小極大歸一化,消除數據量綱和數量級影響。

2) 設輸入層神經元有5個輸入節點,輸出層神經元有2個輸出節點,隱含層神經元個數預設為20個。這樣,輸入層到隱含層間的連接權值w就為20×5的矩陣,權值矩陣的元素一般是[-1,1]之間的隨機數。

3)記b為隱含層的偏置,對應隱含層的L=20個神經元節點,則偏置向量b的大小為20×1的列向量。偏置向量的元素一般是[0,1]之間的隨機數。

4)記β為隱含層與輸出層的連接權值,β大小為20×2的矩陣,當輸出層的神經元節點m=2時,訓練得到的是一個雙輸出的極限學習機模型。通過初始的隱含層權值w和偏置b,即可求解出輸出層的權值矩陣β,從而得到訓練好的極限學習機模型。

5) 通過激活函數G(·)進行特征映射,可以得到極限學習機隱含層的輸出值為:

(2)

6) 設T為期望的輸出值,從而只需要求得合適的β使得誤差函數值最小或者接近于0:

Hβ=T

(3)

對上面的等式求逆,解出輸出層的權值β,如下式所示:

β=H+T

(4)

這樣就可得到一個訓練好的極限學習機模型。

在確定了極限學習機的隱含層權值w和偏置b,以及輸出層的權值矩陣β后,代入需要預測的樣本特征x,即可獲得相應的ELM預測輸出值。

2.3 斐波那契法尋優隱含層節點數

極限學習機參數的隨機初始化使得ELM具有較好的泛化性,但也要求極限學習機模型增加較多的節點數目來實現準確的訓練。在大樣本情況下,過多的節點會消耗計算資源,也可能造成過擬合。

為了解決這個問題,將如何確定ELM的隱含層節點數轉換為了一維約束優化問題,并采用斐波那契搜索法進行求解,稱為斐波那契極限學習機(FELM)。給定一個隱含層節點數的區間,斐波那契搜索算法就可以同時改變搜索區間的2個端點,不斷縮短包含極值點的搜索區間,直到一定程度,這是一種雙向收縮法,其主要優點是可以得到逼近全局的最優解,同時算法的時間復雜度低,修改參數簡單。

執行的步驟如下:

1) 選取初始隱含層節點數,下限a1=5,上限b1=50,由此確定搜索區間為[5,50],由于隱含層個數均為整數,因此給出搜索步長ε=1,求出搜索次數n,使得:

(5)

其中:Fn為斐波拉契數列,F0=F1=1,并設定搜索精度δ=0.05。

2) 當k=1時,計算最初2個搜索點:

(6)

3) 如果f(p1)>f(q1)時,

(7)

4) 如果f(p1)

(8)

5) 當進行至k=n-2時,就無法借比較函數值f(p1)和f(q1)的大小來確定最終區間。為此得到:

pn=pn-1,qn=pn-1+δ

(9)

其中:δ為搜索精度,在pn和qn這兩點中,以兩者的平均值為近似極小值點,相應的函數值為近似極小值。本文通過上述斐波那契對稱搜索得到隱含層個數為24,求解次數僅為9步,優化速度非常快,適合于工業數據建模,因此替換3.2節的第三步的隱含層個數,后面的偏置和權值個數也做相應改變。

3 平衡優化器

平衡優化器(equilibrium optimizer,EO)是于2020年提出的一種源自動態質量平衡方程的新型優化算法,具有尋優能力強,收斂速度快的特點[15]。

概括來說,控制容積內的濃度變化率主要由容量流率和平衡狀態下的濃度增量,以及控制容積內部的質量生成速率這三者決定,因此有如下公式描述:

(10)

其中:Ceq-C是平衡狀態下的濃度增量;Q是容積的容量流率;G是質量生成速率;V是容積;C是容積內濃度。

求解上述微分方程,可得平衡優化器的尋優公式如下:

(11)

(12)

其中:F是指數項;C0是容積內初始濃度;λ是流動速率。

3.1 EO優化步驟

根據式(11)平衡優化器開始迭代尋優,主要步驟描述如下:

步驟1參數初始化。

在初始化步驟中,主要工作包括EO算法的參數設置和生成粒子初始濃度。

1) 參數初始化:初始化迭代計數器Iter,最大迭代次數Max_iter,解空間的邊界限制,大小為N的種群和問題的維數M。此外,還應給出a和GP,用于更新步驟4中的全局搜索權重和生成概率。

2) 粒子初始化:生成N個粒子,每個粒子的維數為M。每個粒子內部的維數隨機初始化為式(13)。

(13)

其中:Cmin、Cmax分別表示搜索空間的最大和最小界;ri代表在區間[0,1]內隨機生成的一個數。

步驟2評估每個候選粒子的適應度值。

將當前迭代中每個粒子的適應度值與前一次迭代的適應度值進行比較,當達到更好地適應度值時將被保存。

步驟3構建平衡池(Ceq)。

如式(14)和式(15)所示,EO算法選取最優的4個粒子Ceq_1、Ceq_2、Ceq_3、Ceq_4,最優的4個粒子對均衡池的平均值Ceq_ave為平衡的候選解。

Ceq_pool={Ceq_1,Ceq_2,Ceq_3,Ceq_4,Ceq_ave}

(14)

(15)

步驟4更新指數項(F)和生成速率(G)。

指數項F表示為式(16):

F=a·sign(r-0.5)(e-λt-1)

(16)

(17)

其中:a為全局搜索的權重;Iter和Max_iter分別表示當前迭代次數和最大迭代次數。

生成速率G如式(18)所示。

G=G0F

(18)

G0=GCP(Ceq-λC)

(19)

(20)

其中:r1、r2為區間[0,1]內均勻分布的隨機數;GCP為生成速率控制參數,表示生成項對更新過程貢獻的可能性,由式(20)計算;GP為生成概率,表示生成項對更新過程的貢獻概率。

步驟5更新每個粒子的濃度(Cnew)。

各粒子濃度的更新規則描述如下:

(21)

其中:C0為控制體積內濃度;λ為停留時間的倒數;V為單位體積。

步驟6執行步驟3,直到迭代Iter等于Max_iter。

3.2 中心EO算法

該算法中找最好的4個解Ceq_1、Ceq_2、Ceq_3、Ceq_4作為候選集,這就失去了種群的多樣性,容易陷入局部極小。本文采用KMedoids對其改進,采用KMedoids聚類,將平衡池中的全部狀態聚成4個類,這4個類的中心也是平衡池中的4個狀態,就用這4個聚類中心作為候選解Ceq_1、Ceq_2、Ceq_3、Ceq_4,平均候選解Ceq_ave仍然是這4個候選解的平均值。

KMedoids算法比較簡單,首先任意選擇初始代表種群,通過用非代表種群替換代表種群來提高聚類質量。代價函數用于評價聚類質量,其計算公式組如下:

(22)

(23)

Cost=Enew-Eold

(24)

其中:k是聚類個數;p是一個樣本;oi是簇的中心樣本;Enew和Eold分別是新、舊中心離差平方和;Cost是新中心替換舊中心的總代價。

如果當前的代表種群被非代表種群所代替,代價函數就計算絕對誤差值的差,交換的總代價是所有非代表種群所產生的代價之和。如果總代價是負的,說明中心替換后簇內成員更緊密了,實際的絕對誤差將會減小,則可以取代,反之,則本次迭代沒有變化。

4 KEO-FELM表冷器預測控制模型

對空調表冷系統,由主表冷閥門、副表冷閥門和電動三通閥以及室外溫濕度等環境變量組成輸入參數,由斐波拉契ELM模型(FELM)可得到系統d步延遲的預測輸出yp(k+d)。為了進行反饋校正,本文采用下式作為預測偏差:

e(k+d)=y(k+d)-yp(k+d)

(25)

中心平衡優化器(KEO)對主表冷閥門、副表冷閥門和電動三通閥的開度進行優化的適應度函數定義為:

(26)

約束條件:

0

其中:ys(k+d)是設定輸出;yp(k+d)是預測輸出;n是輸出變量的個數;d是預測的步數。對于空調表冷系統,預測控制如圖3所示。

圖3 KEO優化的表冷系統FELM預測控制

將一組主表冷閥門、副表冷閥門和電動三通閥的開度作為KEO要尋優的粒子種群,其適應度取式(26)的函數及其約束條件,在有限的預測步數內,以目標函數最小為評價標準搜索最優的控制向量u(k)=[u1(k),u2(k),u3(k)]。實施步驟如圖4所示。

圖4 預測控制模型的流程框圖

5 仿真研究

5.1 數據采集

表冷器溫濕度預測控制模型的節能效果需要進行評估。采集夏季除濕工況的數據作為評估數據集,各種工況下空調機組的運行模式由PLC程序完成,并由傳感器采集和工業以太網傳輸表冷器開度變化趨勢。數據采集結果表明,卷煙廠卷包車間內所有數據監測點的溫度均在21.5~28.7 ℃范圍內,相對濕度均在61.2%~68.0%內,數據采集精度均達到了工藝要求。

針對本文的空調表冷系統,假設采樣周期為5 s,系統輸出對輸入的延遲d=1,選取3個閥門開度u1(k)、u2(k)和u3(k)以及室外溫濕度o1(k)、o2(k)作為輸入,室內溫濕度y1(k+1)、y2(k+1)的序列作為輸出,按照上述運行環境,采集某卷煙廠卷包車間實際運行數據,所有數據均為數值類型。得到1 000組數據,將其中70%的數據用于訓練FELM預測模型,剩余30%的數據用來測試該模型。

5.2 預測模型性能實驗

本文主要是控制表冷器的3個水閥,具體控制思路是:由于室外相對濕度隨季節的變化,室內平均濕度用于控制主表冷閥開度和電動三通閥開度。同樣,由于室外氣溫的季節性變化,室內平均溫度用來控制副表冷閥開度和電動三通閥開度,確保副表冷閥進水溫度滿足降溫要求。由此可以確定預測控制系統的輸入量為:主表冷閥開度、副表冷閥開度和電動三通閥開度,輸出量為室內平均濕度和室內平均溫度。

采用平均百分比誤差(MAPE)和均方根誤差(RMSE)指標評估預測模型精度。

(27)

(28)

其中:yip(k+d)為第i組樣本的預測值;yi(k+d)為第i組樣本的真實值;n為樣本數量。

ELM模型輸入層節點數為5,輸出層節點數為2,隱含層節點數由Fibonacci搜索得到,相應的偏置值個數由隱含層節點數決定。Fibonacci算法的搜索步長ε=1,搜索精度δ=0.05。

除了和基礎ELM比較,本文算法還與文獻[7]使用回聲狀態網絡ESNs建立的預測模型進行對比,由于文獻[7]的模型性能已經優于BP網絡和LS-SVM網絡,因此本文就沒有再和這兩者進行對比了。按照文獻[7]的回聲狀態網絡參數設置,儲備池規模N=200,稀疏度XD=0.03,ρ(W)=0.8,權值Win、W和Wfd隨機產生。

針對測試數據,對比3種預測模型的性能,3種模型輸出的室內溫度和室內相對濕度預測結果如圖5和圖6所示。

從圖5和圖6可以看出,FELM模型預測誤差普遍比ELM和ESNs減少,這表明模型具有較高的預測精度。

3種模型的平均百分比誤差(MAPE)和均方根誤差(RMSE)對比如表1所示。

同樣,由表1這2個指標分析可知,FELM模型的預測效果較ELM和ESNs更好。

圖5 室內溫度FELM模型預測結果

圖6 室內濕度 ELM模型預測結果

表1 模型性能

5.3 控制性能實驗

在仿真控制中,提出了基于KEO優化的表冷系統FELM預測控制的設定。

平衡優化算法的種群數目為30個,最大迭代次數為500次,閥門開度優化的上限值為100%,下限值為0%,全局搜索的權重a為2。KMedoids聚類算法的中心數設定為4,對應于4個候選解。

本文算法與常規基于極限學習機預測控制算法(EO-ELM-MPC)進行了比較,EO-ELM-MPC的參數設置與KEO-FELM-MPC中原生部分是相同的,本文算法也和文獻[7]中PSO-ESNs模型進行了對比。為了測試預測控制模型的階躍響應能力,設定的室內溫度、室內相對濕度分別為24 ℃和64%。預測控制模型的階躍響應如圖7所示。表2和表3展示了溫濕度控制性能指標。

圖7 溫濕度階躍響應曲線

表2 室內溫度控制性能指標

表3 室內濕度控制性能指標

由圖7以及表2和表3可以看出:KEO-FELM-MPC算法相較EO-ELM-MPC和PSO-ESNs-MPC算法在系統的超調量上有明顯降低,具有更好的穩定性,但在上升時間上會有一定的提高,說明對于原有算法的改進增加了計算量,在瞬態性能方面有所降低,但是最終的調節時間并沒有多少延遲,因此該預測控制系統總體性能得到了提高。

為了測試預測控制模型的跟蹤能力,設置的系統參考軌跡為:

(29)

(30)

其中:t為當前采樣時刻;T為總的測試時段。控制系統的跟蹤能力測試如圖8所示。

從圖8可知,對于空調表冷系統,當參考軌跡發生變化時,基于KEO-FELM預測控制算法(KEO-FELM-MPC)能夠快速地優化出主表冷閥門、副表冷閥門和電動三通閥的開度,使空調系統能夠及時跟蹤室內溫度和相對濕度的設定值,而常規EO-ELM-MPC和PSO-ESNs-MPC算法雖然最終也能跟蹤設定值的變化,但是每次跟蹤有躍變的設定值時,產生的超調量較大,調節時間更長,造成控制效果相對較差。因此本文中提出的KEO-FELM-MPC模型的跟蹤性能優于EO-ELM-MPC和PSO-ESNs-MPC模型。

5.4 節能效果對比

通過KEO優化的FELM預測控制算法獲得控制量,即主表冷閥開度、副表冷閥開度和電動三通閥開度,閥門開度變化曲線如圖9所示。

圖8 表冷系統跟蹤設定值曲線

圖9 閥門開度變化曲線

圖9顯示了3個閥門開度變化,圖中KEO-FELM-MPC算法得到的閥門開度均小于EO-ELM-MPC和PSO-ESNs-MPC方法,這說明只要三者協同工作,采用很小的控制量就可以達到較好的控制效果。為了進一步衡量節能效果,本文將閥門控制開度進行了對比,如表4所示。

表4 閥門開度

由表4可知,按照這樣的控制策略調節表冷器閥門,表冷器可以快速準確地達到控制目標,而且閥門開度整體變小,從而減少了冷水供應量,實現了空調系統的節能降耗目的。

6 結論

采用斐波拉契搜索法尋優到了極限學習機隱含層的節點個數,提出了斐波拉契極限學習機(FELM),提高了預測精度;采用KMedoids聚類增強了平衡優化算法的種群粒子的多樣性,提出了中心平衡優化算法(KEO),提高了全局尋優的能力;針對多變量非線性表冷系統,將這2種改進的方法用于表冷系統控制的預測模型和滾動優化中,提出一種基于KEO優化的FELM預測控制算法,根據系統的輸入輸出數據建立FELM預測模型,結合原本的反饋誤差對預測模型進行校正,定義了優化算法的適應度函數,利用KEO算法滾動優化獲得最優控制量,即表冷系統3個閥門的開度。仿真實驗表明,該預測控制模型具有更高的穩定性和跟蹤性,而且具有更好的節能效果。

下一步的工作可以采用帶有時間記憶功能的深度學習預測算法,比如LSTM或GRU建立預測模型,研究如何使這些高性能的深度網絡能快速地建立預測模型,還可以進一步研究控制性能的多目標滾動優化問題。

猜你喜歡
閥門優化模型
一半模型
美嘉諾閥門(大連)有限公司
流程工業(2022年3期)2022-06-23 09:41:08
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
裝配式玻璃鋼閥門井的研發及應用
煤氣與熱力(2021年3期)2021-06-09 06:16:18
一道優化題的幾何解法
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲国产成人久久77| 伊人婷婷色香五月综合缴缴情| 亚洲福利网址| 亚洲精品麻豆| 亚洲av色吊丝无码| 香蕉色综合| 免费观看无遮挡www的小视频| 亚洲精品无码av中文字幕| 91美女视频在线观看| 五月婷婷精品| 专干老肥熟女视频网站| 亚洲天堂视频在线播放| 999福利激情视频| 亚洲精品自产拍在线观看APP| 婷婷色婷婷| 欧美狠狠干| 激情成人综合网| 久久人午夜亚洲精品无码区| 日本影院一区| a毛片在线| 久久精品国产在热久久2019| 国产高清免费午夜在线视频| 国产69囗曝护士吞精在线视频| 欧美第一页在线| 18禁不卡免费网站| 国产后式a一视频| 亚洲精品动漫| 午夜毛片免费观看视频 | 亚洲国产欧美国产综合久久 | 亚洲69视频| 久久一本精品久久久ー99| 国产精品三区四区| 久草网视频在线| 久久亚洲国产视频| 日本高清免费不卡视频| 亚洲人成网线在线播放va| 亚洲精品片911| 国产亚洲视频播放9000| 久热re国产手机在线观看| 欧美性天天| 美女扒开下面流白浆在线试听 | 亚洲欧洲日产国码无码av喷潮| 91精品国产自产在线老师啪l| 国产精品lululu在线观看| 丁香婷婷激情网| 自拍亚洲欧美精品| 色网在线视频| 四虎免费视频网站| 亚洲水蜜桃久久综合网站| 欧美无专区| 欧美色伊人| 99久久精品免费观看国产| 欧美精品一二三区| 亚洲AV无码久久精品色欲| 九月婷婷亚洲综合在线| 91香蕉国产亚洲一二三区| 亚洲综合天堂网| 婷婷丁香色| 综合社区亚洲熟妇p| 五月天久久综合| 天堂成人在线视频| 丝袜亚洲综合| 亚洲一区波多野结衣二区三区| 伊在人亚洲香蕉精品播放| 亚洲免费福利视频| 91精品国产91久无码网站| 亚洲浓毛av| 97国产精品视频人人做人人爱| www亚洲天堂| 国产成人综合久久精品下载| 91蝌蚪视频在线观看| 91视频国产高清| 色成人综合| 亚洲日本一本dvd高清| AV色爱天堂网| 日韩一级二级三级| 九九香蕉视频| 91福利免费视频| 伊人无码视屏| 欧美精品黑人粗大| 欧美国产综合色视频| 久久免费视频播放|