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

高爐煉鐵過程多元鐵水質(zhì)量非線性子空間建模及應(yīng)用

2016-12-17 08:23:46宋賀達(dá)周平王宏柴天佑
自動(dòng)化學(xué)報(bào) 2016年11期
關(guān)鍵詞:方法質(zhì)量模型

宋賀達(dá) 周平 王宏,2 柴天佑

高爐煉鐵過程多元鐵水質(zhì)量非線性子空間建模及應(yīng)用

宋賀達(dá)1周平1王宏1,2柴天佑1

高爐煉鐵是一個(gè)物理化學(xué)反應(yīng)復(fù)雜、多相多場耦合的大滯后、非線性動(dòng)態(tài)系統(tǒng),其關(guān)鍵工藝指標(biāo)—鐵水質(zhì)量參數(shù)的檢測、建模和控制一直是冶金工程和自動(dòng)控制領(lǐng)域的難題.本文提出一種面向控制的數(shù)據(jù)驅(qū)動(dòng)高爐煉鐵多元鐵水質(zhì)量非線性子空間建模方法.首先,為了提高建模效率和降低計(jì)算復(fù)雜度,采用數(shù)據(jù)驅(qū)動(dòng)典型相關(guān)性分析與相關(guān)性分析相結(jié)合的方法提取與鐵水質(zhì)量相關(guān)性最強(qiáng)的關(guān)鍵可控變量作為建模的輸入變量;同時(shí),為了更好地反映高爐非線性動(dòng)態(tài)特性,將相關(guān)輸入輸出變量的時(shí)序和時(shí)滯關(guān)系在建模過程進(jìn)行考慮;最后,采用基于最小二乘支持向量機(jī)(Least square support vector machine, LS-SVM)的非線性Hammerstein系統(tǒng)子空間辨識(shí)方法建立數(shù)據(jù)驅(qū)動(dòng)的多元鐵水質(zhì)量非線性狀態(tài)空間模型.同時(shí),將核函數(shù)表示的模型非線性特性用多項(xiàng)式函數(shù)擬合,在僅損失很小模型精度的前提下大大降低模型的計(jì)算復(fù)雜度.基于實(shí)際數(shù)據(jù)的工業(yè)試驗(yàn)驗(yàn)證了所提建模方法的準(zhǔn)確性、有效性和先進(jìn)性.

高爐煉鐵,子空間辨識(shí),Hammerstein模型,典型相關(guān)性分析,非線性系統(tǒng)建模,最小二乘支持向量機(jī),多元鐵水質(zhì)量

鋼鐵是國民經(jīng)濟(jì)的命脈,高爐煉鐵則是鋼鐵生產(chǎn)中最為重要的環(huán)節(jié)之一.如圖1所示,整個(gè)高爐煉鐵系統(tǒng)分為高爐本體、給料系統(tǒng)、熱風(fēng)系統(tǒng)、煤粉噴吹系統(tǒng)、高爐煤氣處理系統(tǒng)以及出鐵系統(tǒng)等幾個(gè)子系統(tǒng),其中高爐本體又可分為爐喉、爐身、爐腰、爐腹和爐缸等幾個(gè)部分.高爐煉鐵時(shí),鐵礦石、焦炭、溶劑按一定比例及布料制度逐層從高爐頂部裝載到爐喉中,在其下部將預(yù)熱的空氣、氧氣和煤粉鼓入爐缸中,空氣、氧氣、煤粉和焦炭在高溫作用下發(fā)生一系列復(fù)雜物理化學(xué)反應(yīng),生成大量的高溫還原性氣體,這些還原性氣體不斷向上運(yùn)動(dòng)將鐵從鐵礦石中還原出來,而上行的氣體最終變?yōu)楦郀t煤氣從爐頂回收.爐料則隨著爐缸中焦炭的不斷燃燒和鐵水的不斷滴落逐漸向下運(yùn)動(dòng),在下降過程中,爐料經(jīng)過加熱、還原、熔化等一系列復(fù)雜的物理化學(xué)變化,最終生成液態(tài)的生鐵和爐渣從出鐵口排出[1?2].

圖1 高爐煉鐵過程示意圖Fig.1 diagram of a typical BF ironmaking process

高爐煉鐵的目標(biāo)是實(shí)現(xiàn)冶煉過程的高產(chǎn)量與低能耗,為了實(shí)現(xiàn)這一目標(biāo),就應(yīng)對(duì)高爐內(nèi)部狀態(tài)進(jìn)行實(shí)時(shí)監(jiān)測與有效控制.然而高爐內(nèi)部冶煉環(huán)境極其嚴(yán)酷,反應(yīng)最劇烈的區(qū)域溫度高達(dá)2000多度,壓強(qiáng)高達(dá)標(biāo)準(zhǔn)大氣壓的4倍左右,且伴隨著固、液、氣多相共存的狀態(tài),使高爐內(nèi)部狀態(tài)的實(shí)時(shí)監(jiān)測難以實(shí)現(xiàn),從而難以對(duì)高爐進(jìn)行優(yōu)化控制.目前,被廣泛用來間接反映高爐內(nèi)部狀態(tài)的指標(biāo)為鐵水質(zhì)量參數(shù),綜合性的鐵水質(zhì)量指標(biāo)通常采用硅含量([Si])、磷含量([P])、硫含量([S])和鐵水溫度(Molten iron temperature,MIT)來衡量.硅含量是評(píng)定高爐爐況穩(wěn)定性和生鐵質(zhì)量的重要指標(biāo),也是表征高爐熱狀態(tài)及其變化的標(biāo)志之一,硅含量的變化可由多種高爐內(nèi)部狀態(tài)變化引起,比如,熔化區(qū)域的垂直位移、原料的變更和熱量流動(dòng)等[3].硫和磷在鋼材中均是有害元素,在高爐冶煉過程中應(yīng)嚴(yán)格控制鐵水中的硫含量和磷含量.硫含量過高說明渣鐵之間的脫硫反應(yīng)不夠徹底.鐵水溫度是物理熱的表征,是影響高爐穩(wěn)定順行及節(jié)能降耗的重要指標(biāo).采用這四種鐵水質(zhì)量參數(shù)作為高爐內(nèi)部狀態(tài)的評(píng)判指標(biāo)可以較全面地了解高爐內(nèi)部的運(yùn)行狀態(tài),為高爐的控制操作提供指導(dǎo).因此,要實(shí)現(xiàn)面向高爐鐵水質(zhì)量和穩(wěn)定順行的控制與優(yōu)化,首先要建立鐵水質(zhì)量的動(dòng)態(tài)模型.

高爐煉鐵系統(tǒng)是一個(gè)物理化學(xué)反應(yīng)復(fù)雜、多相多場耦合的非線性、大滯后、動(dòng)態(tài)時(shí)變系統(tǒng),因此,鐵水質(zhì)量的預(yù)測模型應(yīng)該是一個(gè)非線性的動(dòng)態(tài)模型.目前,針對(duì)鐵水質(zhì)量建模的研究已有很多,包括基于機(jī)理分析的數(shù)學(xué)模型[4]、基于規(guī)則的推理模型[5]以及基于數(shù)據(jù)的統(tǒng)計(jì)模型[6].近些年來,隨著計(jì)算機(jī)的快速發(fā)展,數(shù)據(jù)驅(qū)動(dòng)的建模方法已被廣泛應(yīng)用于解決鐵水質(zhì)量的建模問題.數(shù)據(jù)驅(qū)動(dòng)是一種黑箱建模方法,它的主要思想是利用一些數(shù)學(xué)工具和智能算法基于過程數(shù)據(jù)直接近似得到過程的輸入輸出關(guān)系而不需要任何先驗(yàn)知識(shí)[7?9].已有的鐵水質(zhì)量數(shù)據(jù)驅(qū)動(dòng)建模方法有偏最小二乘法[10]、人工神經(jīng)網(wǎng)絡(luò)(Artificial neural network,ANN)[11?12]、支持向量回歸機(jī)(Support vector regression,SVR)[13?14]等.但是,這些鐵水質(zhì)量預(yù)測模型大多只是對(duì)單一鐵水質(zhì)量參數(shù)硅含量或鐵水溫度的建模,并不能全面地反映高爐內(nèi)部復(fù)雜的狀態(tài),且部分模型并未考慮高爐系統(tǒng)的非線性與動(dòng)態(tài)特性.另外,這些以鐵水質(zhì)量軟測量為目標(biāo)的輸入輸出智能模型用于解決鐵水質(zhì)量尤其是多元鐵水質(zhì)量控制系統(tǒng)設(shè)計(jì)時(shí),總是存在諸多局限性和不足.雖然,Zeng等[15]采用子空間辨識(shí)方法建立了鐵水硅含量的線性輸入輸出預(yù)測模型,Marutiram等[16]采用線性回歸的方法建立了鐵水硅含量的線性ARMA(Auto regressive and moving average)預(yù)測模型,并且二者均針對(duì)建立的線性模型進(jìn)行了預(yù)測控制,取得了一定的效果.但是,高爐本身是一個(gè)極其復(fù)雜的非線性動(dòng)態(tài)系統(tǒng),用線性模型去近似非線性系統(tǒng),并不能完全準(zhǔn)確地表征鐵水質(zhì)量,因而,基于線性模型的線性控制方法也難于對(duì)多元鐵水質(zhì)量指標(biāo)進(jìn)行有效控制.

針對(duì)上述問題,為了建立一種既適合于控制器設(shè)計(jì),又能充分體現(xiàn)高爐非線性動(dòng)態(tài)特性,同時(shí)又無需任何先驗(yàn)知識(shí)的高爐煉鐵多元鐵水質(zhì)量非線性預(yù)測模型,本文首先采用典型相關(guān)性分析與相關(guān)性分析相結(jié)合的方法選取與鐵水質(zhì)量相關(guān)性最強(qiáng)的高爐本體可控變量作為模型的輸入變量,然后采用基于LS-SVM(Least square support vector machine)的非線性子空間辨識(shí)方法[17]建立數(shù)據(jù)驅(qū)動(dòng)的多元鐵水質(zhì)量非線性狀態(tài)空間動(dòng)態(tài)模型.同時(shí),將核函數(shù)表示的非線性特性用多項(xiàng)式函數(shù)擬合,以降低模型的復(fù)雜度.最后基于實(shí)際工業(yè)現(xiàn)場數(shù)據(jù)進(jìn)行工業(yè)試驗(yàn)和比較分析.狀態(tài)空間模型是現(xiàn)代控制理論的經(jīng)典模型,針對(duì)此模型的控制方法研究已有很多,且理論也很成熟.另外,Hammerstein模型由于其特殊的線性與非線性模塊串聯(lián)結(jié)構(gòu),相較于SVR和ANN等非線性智能模型,更易于設(shè)計(jì)非線性控制器.所以本文建立的多元鐵水質(zhì)量模型可以直接應(yīng)用于多元鐵水質(zhì)量的控制研究中,是一種面向控制的模型.

1 建模策略

從提高產(chǎn)品質(zhì)量、節(jié)約能源的角度而言,高爐系統(tǒng)控制與優(yōu)化的主要對(duì)象是鐵水質(zhì)量參數(shù),尤其是鐵水[Si]和鐵水溫度,它們也是衡量高爐內(nèi)熱狀態(tài)和穩(wěn)定順行的重要參數(shù),其過高和過低對(duì)于燃料消耗和成本有較大的影響.為此,對(duì)鐵水質(zhì)量參數(shù)進(jìn)行建模意義重大,這也是實(shí)現(xiàn)鐵水質(zhì)量控制與優(yōu)化的關(guān)鍵.針對(duì)前述已有鐵水質(zhì)量建模存在的諸多問題,提出圖2所示的多元鐵水質(zhì)量非線性動(dòng)態(tài)建模策略:

圖2 多元鐵水質(zhì)量參數(shù)非線性子空間建模策略Fig.2 Strategy diagram of multivariable nonlinear dynamic modeling for molten iron quality

1)高爐本體參數(shù)較多,且變量之間存在著較強(qiáng)的相關(guān)性,若所有的變量均被用來建立鐵水質(zhì)量模型,勢必會(huì)增大建模過程的計(jì)算復(fù)雜度,影響模型的準(zhǔn)確性和有效性.已有的鐵水質(zhì)量建模研究中,輔助變量的選取方法有機(jī)理分析、主成分分析(Principle component analysis,PCA)和相關(guān)性分析等[11?12,18?19],但是這些選取方法大多是針對(duì)單一鐵水質(zhì)量,并未考慮對(duì)多個(gè)鐵水質(zhì)量指標(biāo)的綜合影響,且已有研究中選取的變量多是針對(duì)建模問題的,在控制問題中有些變量并不適用.所以,采用典型相關(guān)性分析方法選擇與鐵水質(zhì)量指標(biāo)相關(guān)性最大的幾個(gè)可控變量,并用相關(guān)性分析方法對(duì)這些可控變量做進(jìn)一步篩選,去除冗余變量.

2)子空間辨識(shí)建模是一種由輸入輸出數(shù)據(jù)辨識(shí)系統(tǒng)狀態(tài)空間模型的方法,它有三種基本算法,分別為N4SID(Numerical algorithms for subspace state space system identification)[20]、MOESP (Multi-variable output error state space)[21]、CVA (Canonical variate analysis)[22].子空間辨識(shí)可以在系統(tǒng)狀態(tài)完全未知的情況下,僅由輸入輸出數(shù)據(jù)將系統(tǒng)的狀態(tài)估計(jì)出來,從而構(gòu)造系統(tǒng)的狀態(tài)空間模型.本文采用基于LS-SVM的非線性Hammerstein系統(tǒng)子空間建模方法[17]建立多元鐵水質(zhì)量參數(shù)的狀態(tài)空間模型.該方法在傳統(tǒng)線性子空間建模算法N4SID基礎(chǔ)上,引入LS-SVM去解決無需任何先驗(yàn)知識(shí)的靜態(tài)非線性環(huán)節(jié)的辨識(shí)問題,最終可建立一個(gè)能夠較為準(zhǔn)確刻畫煉鐵過程動(dòng)態(tài)特性和非線性特性的狀態(tài)空間模型.該模型不但能準(zhǔn)確地對(duì)多元鐵水質(zhì)量進(jìn)行在線動(dòng)態(tài)估計(jì),并且可用于解決面向多元鐵水質(zhì)量的控制和優(yōu)化問題.

3)由于基于LS-SVM非線性Hammerstein子空間模型的非線性特性是用核函數(shù)表示,計(jì)算效率較低.因此,采用多項(xiàng)式擬合方法[23]擬合模型的非線性特性,即將核函數(shù)表示的非線性部分替換為計(jì)算較為簡單的多項(xiàng)式表達(dá)式,以降低模型的計(jì)算復(fù)雜度.

2 建模算法

2.1 基于典型相關(guān)分析和相關(guān)性分析的建模輔助變量選擇

典型相關(guān)分析的基本思想與PCA類似,首先在兩組變量(2)所示,(y1,y2,···,yq)中分別找出變量的一個(gè)線性組合,即

式中,α(i),β(i)分別為線性組合的系數(shù)向量,使得兩組線性組合之間具有最大的相關(guān)系數(shù).然后選取相關(guān)系數(shù)僅次于第一對(duì)線性組合并且與第一對(duì)線性組合不相關(guān)的第二對(duì)線性組合,如此繼續(xù)下去,直到兩組變量之間的相關(guān)性被提取完畢為止.Ui,Vi的配對(duì)稱為典型變量,它們之間的相關(guān)系數(shù)稱為典型相關(guān)系數(shù).

在進(jìn)行輔助變量選擇時(shí),首先要對(duì)各對(duì)典型變量的典型相關(guān)系數(shù)進(jìn)行顯著性檢驗(yàn),若某一對(duì)典型變量的相關(guān)程度不顯著,說明這對(duì)變量不具有代表性,舍棄這一對(duì)典型變量.因?yàn)閁i中系數(shù)絕對(duì)值較大的n個(gè)變量對(duì)Ui起決定性作用,所以若Ui,Vi的相關(guān)程度顯著且又具有很大的相關(guān)系數(shù),則這n個(gè)變量便與Vi,即中q個(gè)變量的線性組合具有很大的相關(guān)性.因此,選取輔助變量時(shí),先找出相關(guān)程度顯著且典型相關(guān)系數(shù)較高的幾對(duì)典型變量,然后選取這幾個(gè)Ui中線性組合系數(shù)絕對(duì)值較大的幾個(gè)變量作為候選輔助變量.

在高爐本體變量中,各個(gè)變量均具有一定的相關(guān)性,若兩個(gè)變量相關(guān)性達(dá)到80%以上,則可近似認(rèn)為二者呈線性關(guān)系,不必將二者同時(shí)選為輔助變量.因此,在確定輔助變量時(shí),要用相關(guān)性分析的方法對(duì)候選輔助變量進(jìn)一步縮減,縮減的原則是保留可控變量,舍棄不可控變量.

2.2 基于LS-SVM的非線性Hammerstein系統(tǒng)子空間建模

Hammerstein模型是一類非線性系統(tǒng),由靜態(tài)非線性環(huán)節(jié)和動(dòng)態(tài)線性環(huán)節(jié)串聯(lián)而成.這種模型能較好地反映過程特征,可以描述一大類非線性系統(tǒng),如無線發(fā)射器[24]、聚丙烯牌號(hào)切換系統(tǒng)[25]、燃料電池[26]、生物系統(tǒng)[27]等.在狀態(tài)空間模型輸入端加入靜態(tài)非線性函數(shù)即構(gòu)成Hammerstein模型,如式

在文獻(xiàn)[17]中,Goethals等以預(yù)測為目的提出了一種針對(duì)如式(2)所示的多輸入多輸出Hammerstein系統(tǒng)的系統(tǒng)辨識(shí)方法.通過將N4SID算法中的斜投影計(jì)算問題改寫為一系列分量形式的LS-SVMs回歸問題,從而對(duì)N4SID算法進(jìn)行非線性擴(kuò)展.Goethals等指出Hammerstein系統(tǒng)中的線性模型和式(2)中的靜態(tài)非線性特性可以通過對(duì)LS-SVMs回歸問題得到的矩陣進(jìn)行低秩逼近得到.具體建模過程描述如下.

2.2.1 計(jì)算斜投影

式(4)亦可寫為

式中,s=1,···,li,t=1,···,j.一旦式(4)、(5)中的的估計(jì)值通過最小二乘法確定下來,斜投影可計(jì)算如下:

然而式(5)中包含參數(shù)矩陣Lu與非線性函數(shù)f的相乘項(xiàng),使得優(yōu)化問題非凸,為此引入一組函數(shù)gh,s,k:R→R,

式(10)的形式可以方便地引入核函數(shù)及分量形式(Componentwise)LS-SVM的方法去估計(jì)參數(shù).

式中,ψ:Rn→Rn,δu,δy均為常數(shù)向量.由于系統(tǒng)模型加入了一個(gè)新的參數(shù)δy,因此式(10)變換為

式中?表示Kronecker積.對(duì)于所有h=1,···,2i, s=1,···,li,均有則原LS-SVM回歸問題可表示為一個(gè)有約束的優(yōu)化問題

通過構(gòu)造式 (12)的拉格朗日函數(shù) (Lagrangian),Ly和δy的估計(jì)值可以從如下對(duì)偶系統(tǒng)求解

2.2.2 系統(tǒng)狀態(tài)估計(jì)

2.2.3 提取系統(tǒng)矩陣和靜態(tài)非線性函數(shù)

系統(tǒng)矩陣和靜態(tài)非線性可由式(17)估計(jì)得出

如同上節(jié),這個(gè)最小二乘問題依然可以變?yōu)橐粋€(gè)LS-SVM回歸問題,定義

式中,E為式(17)的殘差.則LS-SVM回歸的原始問題可以寫為

式中,γBD為一個(gè)正則化常數(shù),與第2.2.1節(jié)中的γ相區(qū)別.

式(18)所示ΘAC中的系統(tǒng)矩陣A,C可由如下所示的對(duì)偶問題求出:

對(duì)ΘBD中的系統(tǒng)矩陣B,D的第k列和非線性特性的估計(jì)可通過奇異值分解(SVD)得到.

2.3 多項(xiàng)式擬合模型靜態(tài)非線性特性

然后,將多項(xiàng)式曲線擬合問題轉(zhuǎn)化為如下最小二乘問題求解

利用多元函數(shù)求極值方法,分別對(duì)ai求偏導(dǎo),可得m+1個(gè)方程,寫成矩陣形式如式(22)所示.

求解如式(22)所示的關(guān)于a0,a1,···,am的線性方程組,即可獲得非線性特性的多項(xiàng)式函數(shù)表達(dá)式.從而,Hammerstein模型中用核函數(shù)表示的靜態(tài)非線性特性可由多項(xiàng)式函數(shù)取代,降低了模型的復(fù)雜度.

2.4 建模算法實(shí)現(xiàn)步驟

基于前述討論,最終的輸入變量選擇及非線性子空間辨識(shí)建模過程可以總結(jié)如下:

1)根據(jù)式(1)針對(duì)輸入輸出數(shù)據(jù)求典型相關(guān)性系數(shù).

2)進(jìn)行顯著性檢驗(yàn)和相關(guān)性分析,以此確定建模輸入變量.

3)根據(jù)式(15)和式(16)計(jì)算斜投影Oi和Oi+1.

4)通過對(duì)斜投影Oi,Oi+1進(jìn)行SVD分解估計(jì)出系統(tǒng)的狀態(tài).

5)根據(jù)式(19)獲得系統(tǒng)矩陣A,C及ABD,BBD的估計(jì)值.

6)通過對(duì)式(20)進(jìn)行秩為1的近似,獲得系統(tǒng)矩陣B,D的第k列及非線性特性

3 工業(yè)實(shí)驗(yàn)

3.1 模型建立

本文基于柳鋼2#高爐的實(shí)際生產(chǎn)數(shù)據(jù)進(jìn)行工業(yè)試驗(yàn),建立多元鐵水質(zhì)量模型.由生產(chǎn)現(xiàn)場數(shù)據(jù)庫可得18個(gè)高爐主體參數(shù)分別為:送風(fēng)比、熱風(fēng)壓力、頂壓、壓差、頂壓風(fēng)量比、透氣性、阻力系數(shù)、熱風(fēng)溫度、富氧流量、富氧率、設(shè)定噴煤量、鼓風(fēng)濕度、理論燃燒溫度、標(biāo)準(zhǔn)風(fēng)速、實(shí)際風(fēng)速、鼓風(fēng)動(dòng)能、爐腹煤氣量、爐腹煤氣指數(shù).在建立高爐多元鐵水質(zhì)量模型之前,首先要從這18個(gè)高爐主體參數(shù)中確定出對(duì)多元鐵水質(zhì)量影響最大的幾個(gè)輔助變量.應(yīng)用第2.1節(jié)基于典型相關(guān)分析和相關(guān)性分析的方法進(jìn)行建模輸入變量選擇.首先,采用典型相關(guān)分析的方法針對(duì)柳鋼2#高爐2013年10月10日到18日的200組(采樣間隔為1h)高爐本體數(shù)據(jù)與鐵水質(zhì)量數(shù)據(jù)進(jìn)行分析,去除相關(guān)性不顯著的1對(duì)典型變量,得到相關(guān)系數(shù)分別為0.715、0.571、0.531的3對(duì)典型變量,典型相關(guān)分析結(jié)果如表1所示.

表1 典型相關(guān)分析結(jié)果Table 1 The results of canonical correlation analysis

選出表1中對(duì)輸出典型變量影響較大的幾個(gè)高爐主體參數(shù)作為候選輔助變量,即壓差、阻力系數(shù)、熱風(fēng)溫度、富氧流量、富氧率、設(shè)定噴煤量、理論燃燒溫度和爐腹煤氣量.然后通過相關(guān)性分析的方法選出候選輔助變量中相關(guān)性較大的幾組中的可調(diào)控變量,如富氧流量與富氧率的相關(guān)系數(shù)為0.999,選取富氧率;理論燃燒溫度和富氧率的相關(guān)系數(shù)為0.78,舍棄理論燃燒溫度.壓差和阻力系數(shù)均為模型間接計(jì)算值,考慮到壓差和阻力系數(shù)均通過熱風(fēng)壓力計(jì)算得出,且熱風(fēng)壓力可控,所以選取熱風(fēng)壓力替代壓差和阻力系數(shù).爐腹煤氣量由富氧流量、設(shè)定噴煤量等可調(diào)控量計(jì)算得來,為間接可調(diào)控量,所以選為輔助變量.最終確定用于建模的可調(diào)可控輔助變量為熱風(fēng)壓力、熱風(fēng)溫度、富氧率、設(shè)定噴煤量和爐腹煤氣量.另外,為了反映高爐的非線性動(dòng)態(tài)特性以及輸入輸出變量的時(shí)序關(guān)系,將上一采樣時(shí)刻輔助變量測量值以及上一采樣時(shí)刻多元鐵水質(zhì)量值,連同當(dāng)前采樣時(shí)刻輸入輔助變量檢測值共同作為動(dòng)態(tài)模型的綜合輸入.

高爐是一個(gè)強(qiáng)噪聲系統(tǒng),因此在建模前首先要對(duì)工業(yè)現(xiàn)場采集到的實(shí)際生產(chǎn)數(shù)據(jù)進(jìn)行數(shù)據(jù)預(yù)處理.針對(duì)由于高爐爐況不穩(wěn)定和檢測儀器不精確造成的跳變數(shù)據(jù),采用異常值檢測算法剔除高爐生產(chǎn)過程中的噪聲尖峰跳變數(shù)據(jù);然后采用移動(dòng)平均濾波算法減弱訓(xùn)練數(shù)據(jù)中的高斯噪聲干擾.之后基于預(yù)處理后的數(shù)據(jù),采用本文所提方法進(jìn)行多元鐵水質(zhì)量參數(shù)的非線性狀態(tài)空間模型建模,并采用遺傳算法優(yōu)化模型可調(diào)參數(shù)σ,γ,γBD.當(dāng)可調(diào)參數(shù)確定為σ=1,γ=0.67,γBD=1時(shí),得到的模型如式(23)所示:

圖3 多項(xiàng)式擬合非線性Fig.3 Polynomial fitting for nonlinearity

3.2 建模及其估計(jì)效果

圖4為建立的非線性狀態(tài)空間模型在多項(xiàng)式擬合非線性特性前后對(duì)多元鐵水質(zhì)量參數(shù)(即[Si],[P], [S]和MIT)的建模效果比較,圖5為建立的非線性狀態(tài)空間模型在多項(xiàng)式擬合非線性特性前后對(duì)實(shí)際數(shù)據(jù)的多元鐵水質(zhì)量參數(shù)估計(jì)效果比較.可以看出多項(xiàng)式擬合非線性特性的式(23)所示模型及用核函數(shù)表示非線性特性的模型得到的多元鐵水質(zhì)量預(yù)測值的大小和變化趨勢幾乎完全一致,兩種表示形式的非線性子空間辨識(shí)建立的多元鐵水質(zhì)量模型均具有非常高的建模和估計(jì)精度,建模誤差以及估計(jì)誤差均非常小,且變化趨勢與其實(shí)際值非常一致,即模型估計(jì)輸出與其實(shí)際值已基本吻合.另外,由于建立的模型既能對(duì)多元鐵水質(zhì)量輸出進(jìn)行準(zhǔn)確估計(jì),并且還能輸出系統(tǒng)的狀態(tài)估計(jì)曲線,因而可更好地用于各種控制算法的研究和設(shè)計(jì).

圖4 基于非線性子空間辨識(shí)的多元鐵水質(zhì)量模型在多項(xiàng)式擬合非線性特性前后的建模結(jié)果Fig.4 Modeling results of molten iron quality prediction with and without polynomial fitting

圖5 模型在多項(xiàng)式擬合非線性特性前后對(duì)多元鐵水質(zhì)量參數(shù)的實(shí)際估計(jì)效果Fig.5 Estimation results of molten iron quality prediction with and without polynomial fitting

為了進(jìn)一步說明多項(xiàng)式擬合非線性特性的模型相比于核函數(shù)表示非線性特性模型的優(yōu)越性,計(jì)算二者在實(shí)際數(shù)據(jù)上預(yù)測的均方根誤差(Root mean square error,RMSE)及計(jì)算時(shí)間進(jìn)行定量分析.在實(shí)際工業(yè)測試數(shù)據(jù)計(jì)算得到的兩種模型的估計(jì)RMSE和計(jì)算時(shí)間如表2所示.可以看出兩種模型的預(yù)測精度相差非常小,核函數(shù)表示非線性特性的模型略好于多項(xiàng)式擬合非線性特性的模型,但是在計(jì)算時(shí)間上多項(xiàng)式擬合非線性特性的模型要明顯優(yōu)于核函數(shù)表示非線性特性的模型,計(jì)算效率提高了十幾倍.因此,綜合以上的分析和評(píng)價(jià),可以看出多項(xiàng)式擬合非線性特性的多元鐵水質(zhì)量非線性Hammerstein模型可以在損失非常小建模精度的前提下,大大降低模型復(fù)雜度,提高實(shí)際應(yīng)用過程中的計(jì)算效率.

表2 兩種模型對(duì)各個(gè)鐵水質(zhì)量指標(biāo)估計(jì)的均方根誤差和計(jì)算時(shí)間比較Table 2 Comparison between two models in RMSE for molten iron quality prediction and computation time

為了進(jìn)一步驗(yàn)證本文方法的優(yōu)越性,將本文建模方法即多項(xiàng)式擬合非線性特性的模型(Nonlinear subspace model,N-SM)與文獻(xiàn)[20]提出的線性子空間辨識(shí)(Linear subspace model,L-SM)方法、文獻(xiàn)[28]提出的多輸出LS-SVR(Multi-output least square vector regression,M-LS-SVR)方法以及常見的RLS-ARMA(Recursive least square-auto regressive and moving average)方法進(jìn)行比較研究.這幾種方法根據(jù)實(shí)際工業(yè)數(shù)據(jù)對(duì)多元鐵水質(zhì)量的估計(jì)效果如圖6所示.可以看出,本文模型得到的各個(gè)鐵水質(zhì)量參數(shù)的估計(jì)值相比其他方法能夠更好地?cái)M合真實(shí)值的變化趨勢,且能保證較高的精度,具有更好的估計(jì)性能.為了進(jìn)一步說明這個(gè)問題,引入誤差自相關(guān)函數(shù)來評(píng)價(jià)不同方法的建模性能.眾所周知,一個(gè)好的模型,它的估計(jì)誤差自相關(guān)函數(shù)曲線應(yīng)該近似為一個(gè)白噪聲序列.所以,畫出不同建模方法得到的模型的估計(jì)誤差自相關(guān)函數(shù)曲線,如圖7所示.可以看出,相對(duì)于L-SM、M-LS-SVR以及RLS-ARMA這幾種算法,本文模型4個(gè)鐵水質(zhì)量參數(shù)估計(jì)誤差自相關(guān)曲線綜合來說更加接近于白噪聲序列.

圖6 不同建模方法對(duì)多元鐵水質(zhì)量的實(shí)際估計(jì)效果對(duì)比Fig.6 Estimation results of molten iron quality prediction with different models

另外,為了進(jìn)一步對(duì)比不同模型的估計(jì)精度,將鐵水質(zhì)量的實(shí)際值與模型估計(jì)值分別作為橫、縱坐標(biāo),畫出上述各個(gè)方法對(duì)各個(gè)指標(biāo)的估計(jì)值與實(shí)際值散點(diǎn)圖,并求出各個(gè)散點(diǎn)圖回歸線的斜率.由圖8可以看出,相較于L-SM、M-LS-SVR以及RLSARMA這幾種建模方法,本文建立的模型得到的各個(gè)鐵水質(zhì)量參數(shù)散點(diǎn)更加集中地分布在對(duì)角線附近,且回歸線的斜率更加接近于1,這進(jìn)一步表明本文模型得到的多元鐵水質(zhì)量估計(jì)值更加接近于各個(gè)指標(biāo)參數(shù)的實(shí)際值.評(píng)價(jià)一個(gè)模型的精度通常還會(huì)選用一些標(biāo)準(zhǔn)的統(tǒng)計(jì)指標(biāo),如均方誤差(MSE)、均方根誤差(RMSE)等.本文選用RMSE作為評(píng)價(jià)指標(biāo)對(duì)4種不同模型的估計(jì)精度進(jìn)行定量評(píng)估.在實(shí)際工業(yè)測試數(shù)據(jù)計(jì)算得到的各個(gè)模型的估計(jì)RMSE如表3所示.可以看出本文模型對(duì)各個(gè)鐵水質(zhì)量指標(biāo)的RMSE數(shù)值均為最小,這進(jìn)一步驗(yàn)證了建立的模型具有最好的估計(jì)精度.

表3 多元鐵水質(zhì)量估計(jì)值均方根誤差比較Table 3 RMSE for molten iron quality prediction

3.3 應(yīng)用于多元鐵水質(zhì)量非線性預(yù)測控制

為了進(jìn)一步驗(yàn)證本文方法建立的多元鐵水質(zhì)量模型的有效性,將建立的Hammerstein狀態(tài)空間模型作為預(yù)測模型應(yīng)用于多元鐵水質(zhì)量的非線性預(yù)測控制.預(yù)測控制的性能指標(biāo)如式(25)所示.

圖7 不同模型估計(jì)誤差自相關(guān)函數(shù)Fig.7 Autocorrelation function of estimating error of different models

選取預(yù)測時(shí)域Np為1,控制時(shí)域Nc為1,控制項(xiàng)權(quán)重R為0.001.多元鐵水質(zhì)量的初始設(shè)定值為[Si]=0.45%,[P]=0.115%,[S]=0.025%, MIT=1500°C,在20、40、60、80時(shí)刻分別更改設(shè)定值,將設(shè)定值分別更改為MIT=1510°C, [S]=0.015%,[P]=0.125%,[Si]=0.50%,在100、120、140、160、180時(shí)刻分別對(duì)熱風(fēng)壓力(u1)、熱風(fēng)溫度(u2)、富氧率(u3)、設(shè)定噴煤量(u4)、爐腹煤氣量(u5)加入較大的脈沖干擾,得到的多元鐵水質(zhì)量響應(yīng)曲線如圖9所示,控制量曲線如圖10所示.可以看出4個(gè)鐵水質(zhì)量參數(shù)均能夠快速地跟隨設(shè)定值的階躍變化,且達(dá)到穩(wěn)態(tài)時(shí)與設(shè)定值的偏差較小,在對(duì)輸入變量加入較大干擾時(shí),鐵水質(zhì)量參數(shù)能夠從較大的波動(dòng)中迅速地回到設(shè)定值,有效地抑制干擾.實(shí)驗(yàn)結(jié)果驗(yàn)證了基于本文模型的非線性預(yù)測控制在多元鐵水質(zhì)量控制上的有效性.

圖8 不同建模方法的估計(jì)值與實(shí)際值散點(diǎn)圖Fig.8 Scatter diagram of estimated and actual value by different models

綜合,可以看出本文建立的模型應(yīng)用于多元鐵水質(zhì)量在線估計(jì)與預(yù)測時(shí),具有較好的估計(jì)性能.另外,本文方法是一種面向控制的非線性系統(tǒng)狀態(tài)空間建模方法,建立的多元鐵水質(zhì)量非線性狀態(tài)空間模型不僅可以用于多元鐵水質(zhì)量參數(shù)的在線估計(jì),為工業(yè)現(xiàn)場的操作人員提供指導(dǎo),并且還可進(jìn)一步應(yīng)用于鐵水質(zhì)量的控制與優(yōu)化.

4 結(jié)論

圖9 多元鐵水質(zhì)量預(yù)測控制結(jié)果Fig.9 Predictive control results of molten iron quality parameters

圖10 多元鐵水質(zhì)量預(yù)測控制控制量曲線Fig.10 Control input curves of predictive control of molten iron quality parameters

高爐煉鐵的最終任務(wù)是生產(chǎn)出符合要求的鐵水.因此,面向[Si]、[P]、[S]和鐵水溫度等多元鐵水質(zhì)量參數(shù)的建模、控制與優(yōu)化歷來是高爐煉鐵生產(chǎn)及其控制工程非常重要的一環(huán).本文針對(duì)現(xiàn)有鐵水質(zhì)量參數(shù)建模的諸多不足,提出一種面向控制的數(shù)據(jù)驅(qū)動(dòng)高爐煉鐵多元鐵水質(zhì)量參數(shù)非線性子空間建模方法.首先,為了提高建模效率和降低計(jì)算復(fù)雜度,采用典型相關(guān)性分析與相關(guān)性分析相結(jié)合的方法從眾多關(guān)聯(lián)變量中提取與鐵水質(zhì)量相關(guān)性最強(qiáng)的5個(gè)關(guān)鍵可控變量作為建模輸入變量;然后基于實(shí)際工業(yè)現(xiàn)場數(shù)據(jù),采用基于LS-SVM的非線性Hammerstein系統(tǒng)子空間辨識(shí)方法,建立數(shù)據(jù)驅(qū)動(dòng)的多元鐵水質(zhì)量非線性狀態(tài)空間動(dòng)態(tài)模型,并將核函數(shù)表示的非線性特性用多項(xiàng)式函數(shù)擬合,以降低模型的復(fù)雜度.工業(yè)試驗(yàn)結(jié)果表明:建立的模型不僅能夠?qū)Χ嘣F水質(zhì)量參數(shù)進(jìn)行準(zhǔn)確在線估計(jì),而且可很好地用于面向鐵水質(zhì)量參數(shù)的控制研究.

1 Gao C H,Jian L,Luo S H.Modeling of the thermal state change of blast furnace hearth with support vector machines.IEEE Transaction on Industrial Electronics,2012, 59(2):1134?1145

2 Kuang S B,Li Z Y,Yan D L,Qi Y H,Yu A B.Numerical study of hot charge operation in ironmaking blast furnace. Minerals Engineering,2014,63:45?56

4 Das S K,Kumari A,Bandopadhay D,Akbar S A,Mondal G K.A mathematical model to characterize effects of liquid hold-up on bosh silicon transport in the dripping zone of a blast furnace.Applied Mathematical Modeling,2011,35(9): 4208?4221

5 Zarandi M H F,Ahmadpour P.Fuzzy agent-based expert system for steel making process.Expert Systems with Applications,2009,36(5):9539?9547

6 Chao Y C,Su C W,Huang H P.The adaptive autoregressive models for the system dynamics and prediction of blast furnace.Chemical Engineering Communications,1986, 44(1?6):309?330

7 Hou Z S,Wang Z.From model-based control to data-driven control:survey,classification and perspective.Information Sciences,2013,235:3?35

8 Hou Z S,Jin S T.Data-driven model-free adaptive control for a class of MIMO nonlinear discrete-time systems.IEEE Transactions on Neural Networks,2011,22(12):2173?2188

9 Zhou P,Lu S W,Chai T Y.Data-driven soft-sensor modeling for product quality estimation using case-based reasoning and fuzzy-similarity rough sets.IEEE Transactions on Automation Science and Engineering,2014,11(4): 992?1003

10 Shi L,Li Z L,Yu T,Li J P.Model of hot metal silicon content in blast furnace based on principal component analysis application and partial least square.Journal of Iron and Steel Research,International,2011,18(10):13?16

11 Zhou P,Yuan M,Wang H,Wang Z,Chai T Y.Multivariable dynamic modeling for molten iron quality using online sequential random vector functional-link networks with self-feedback connections.Information Sciences,2015,325: 237?255

12 Yuan M,Zhou P,Li M L,Li R F,Wang H,Chai T Y. Intelligent multivariable modeling of blast furnace molten iron quality based on dynamic AGA-ANN and PCA.Journal of Iron and Steel Research,International,2015,22(6): 487?495

13 Tang X L,Zhuang L,Jiang C J.Prediction of silicon content in hot metal using support vector regression based on chaos particle swarm optimization.Expert Systems with Applications,2009,36(9):11853?11857

14 Liu Y,Gao Z L.Enhanced just-in-time modelling for online quality prediction in BF ironmaking.Ironmaking and Steelmaking,2015,42(5):321?330

15 Zeng J S,Gao C H,Su H Y.Data-driven predictive control for blast furnace ironmaking process.Computers and Chemical Engineering,2010,34(11):1854?1862

16 Marutiram K,Radhakrishnan V R.Predictive control of blast furnaces.In:Proceedings of the 1991 IEEE International Conference on EC3–Energy,Computer,Communication and Control Systems(TENCON'91).New Delhi,India: IEEE,1991.488?491

17 Goethals I,Pelckmans K,Suykens J A K,De Moor B. Subspace identification of Hammerstein systems using least squares support vector machines.IEEE Transactions on Automatic Control,2005,50(10):1509?1519

19 Phadke M,Wu S M.Identification ofmultiinputmultioutput transfer function and noise model of a blast furnace from closed-loop data.IEEE Transactions on Automatic Control,1974,19(6):944?951

20 Van Overschee P,De Moor B.N4SID:subspace algorithms for the identification of combined deterministic-stochastic systems.Automatic,1994,30(1):75?93

21 Verhaegen M,Dewilde P.Subspace model identification part 1.The output-error state-space model identification class of algorithms.International Journal of Control,1992,56(5): 1187?1210

22 Larimore W E.Canonical variate analysis in identification, filtering,and adaptive control.In:Proceedings of the 29th Conference on Decision Control.Honolulu,HI:IEEE,1990. 596?601

23 Liu Y,Gao Y C,Gao Z L,Wang H Q,Li P.Simple nonlinear predictive control strategy for chemical processes using sparse kernel learning with polynomial form.Industrial& Engineering Chemistry Research,2010,49(17):8209?8218

24 Moon J,Kim B.Enhanced Hammerstein behavioral model for broadband wireless transmitters.IEEE Transactions on Microwave Theory and Techniques,2011,59(4):924?933

25 He De-Feng,Yu Li.Nonlinear predictive control of constrained Hammerstein systems and its research on simulation of polypropylene grade transition.Acta Automatica Sinica,2009,35(12):1558?1563 (何德峰,俞立.約束Hammerstein系統(tǒng)非線性預(yù)測控制及其在聚丙烯牌號(hào)切換中的仿真研究.自動(dòng)化學(xué)報(bào),2009,35(12): 1558?1563)

26 Huo H B,Zhu X J,Hu W Q,Tu H Y,Li J,Yang J.Nonlinear model predictive control of SOFC based on a Hammerstein model.Journal of Power Sources,2008,185(1):338?344

27 Su S W,Wang L,Celler B G,Savkin A V,Guo Y.Identification and control for heart rate regulation during treadmill exercise.IEEE Transactions on Biomedical Engineering,2007,54(7):1238?1246

28 Xu S,An X,Qiao X D,Zhu L J,Li L.Multi-output leastsquares support vector regression machines.Pattern Recognition Letters,2013,34(9):1078?1084

宋賀達(dá) 東北大學(xué)碩士研究生.2013年獲得東北大學(xué)學(xué)士學(xué)位.主要研究方向?yàn)閿?shù)據(jù)驅(qū)動(dòng)建模與控制,機(jī)器學(xué)習(xí)算法.

E-mail:1401124@stu.neu.edu.cn

(SONG He-Da Master student at Northeastern University.He received his bachelor degree from Northeastern University in 2013.His research interest covers data-driven modeling and control and machine learning algorithm.)

周 平 東北大學(xué)副教授.分別于2003年,2006年,2013年獲得東北大學(xué)學(xué)士、碩士和博士學(xué)位.主要研究方向?yàn)楣I(yè)過程運(yùn)行反饋控制,數(shù)據(jù)驅(qū)動(dòng)建模與控制.本文通信作者.

E-mail:zhouping@mail.neu.edu.cn

(ZHOU Ping Associate professor at Northeastern University.He received his bachelor,master and Ph.D.degrees from Northeastern University in 2003,2006 and 2013,respectively.His research interest covers operation feedback control of industrial process,data-driven modeling and control.Corresponding author of this paper.)

王 宏 中組部“千人計(jì)劃”教授,國家特聘專家,IET和InstMC Fellow,長江學(xué)者講座教授,英國曼徹斯特大學(xué)教授.主要研究方向?yàn)榉歉咚闺S機(jī)系統(tǒng)的隨機(jī)分布控制,故障檢測與診斷,非線性控制,基于數(shù)據(jù)的復(fù)雜系統(tǒng)建模.

E-mail:hong.wang@manchester.ac.uk

(WANG Hong IET,InstMC fellow,professor at the Control System Centre,University of Manchester,Manchester,UK since 2002.His research interest covers stochastic distribution control of non-Gaussian stochastic system,fault detection and diagnosis, nonlinear control,and data-based modeling for complex systems.)

柴天佑 中國工程院院士,東北大學(xué)教授,IEEE Fellow,IFAC Fellow.1985年獲得東北大學(xué)博士學(xué)位.主要研究方向?yàn)樽赃m應(yīng)控制,多變量智能解耦控制,流程工業(yè)綜合自動(dòng)化理論、方法與技術(shù).

E-mail:tychai@mail.neu.edu.cn

(CHAI Tian-You Academician of Chinese Academy of Engineering,professor at Northeastern University,IEEE Fellow,IFAC Fellow.He received his Ph.D.degree from Northeastern University in 1985.His research interest covers adaptive control,intelligent decoupling control,and integrated automation theory,method and technology of industrial process.)

Nonlinear Subspace Modeling of Multivariate Molten Iron Quality in Blast Furnace Ironmaking and Its Application

SONG He-Da1ZHOU Ping1WANG Hong1,2CHAI Tian-You1

Blast furnace ironmaking is a nonlinear dynamic process containing complex physical-chemical reaction,multiphase multi-field coupling and large time delay.Measuring,modeling and control of the key process indices of ironmaking process,molten iron quality(MIQ)parameters,have always been treated as a difficult problem in metallurgic engineering and automation field.This paper presents a control oriented data-driven nonlinear subspace modeling method for multivariate prediction of MIQ.First,to improve modeling efficiency,data driven canonical correlation analysis(CCA) and correlation analysis(CA)are combined to pick out the most influential controllable variables from multitudinous factors to serve as the input variables of modeling.Second,to better reflect the nonlinear dynamic characteristics of blast furnace ironmaking process,the time series and time delays of the relevant input and output variables are taken into account.Finally,a data-driven nonlinear state-space model of MIQ is built using least square support vector machine (LS-SVM)based nonlinear subspace identification method for Hammerstein system.With polynomial fitting method,the nonlinear parts expressed by kernel functions in the obtained Hammerstein model are simplified,so as to greatly reduce the computational complexity of the model on the premise of only a small loss of accuracy.Industrial experiments based on real data verifies the accuracy,effectiveness and advancement of the proposed method.

Blast furnace ironmaking,subspace identification,Hammerstein model,canonical correlation analysis (CCA),nonlinear system modeling,least square support vector machine(LS-SVM),multivariate molten iron quality

宋賀達(dá),周平,王宏,柴天佑.高爐煉鐵過程多元鐵水質(zhì)量非線性子空間建模及應(yīng)用.自動(dòng)化學(xué)報(bào),2016,42(11): 1664?1679

Song He-Da,Zhou Ping,Wang Hong,Chai Tian-You.Nonlinear subspace modeling of multivariate molten iron quality in blast furnace ironmaking and its application.Acta Automatica Sinica,2016,42(11):1664?1679

2015-12-04 錄用日期2016-04-18

Manuscript received December 4,2015;accepted April 18,2016

國家自然科學(xué)基金(61473064,61290323,61333007),中央高校基本科研業(yè)務(wù)費(fèi)項(xiàng)目(N130108001),遼寧省教育廳科技項(xiàng)目(L20150186)資助

Supported by National Natural Science Foundation of China (61473064,61290323,61333007),the Fundamental Research Funds for the Central Universities(N130108001),the General Project on Scientific Research for the Education Department of Liaoning Province(L20150186)

本文責(zé)任編委候忠生

Recommended by Associate Editor HOU Zhong-Sheng

1.東北大學(xué)流程工業(yè)綜合自動(dòng)化國家重點(diǎn)實(shí)驗(yàn)室 沈陽110819 中國2.曼切斯特大學(xué)控制系統(tǒng)中心曼徹斯特M60 1QD英國

1.State Key Laboratory of Synthetical Automation for Process Industries,Northeastern University,Shenyang 110819,China 2.Control System Center,Manchester University,Manchester M60 1QD,UK

DOI 10.16383/j.aas.2016.c150819

猜你喜歡
方法質(zhì)量模型
一半模型
“質(zhì)量”知識(shí)鞏固
質(zhì)量守恒定律考什么
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
做夢導(dǎo)致睡眠質(zhì)量差嗎
3D打印中的模型分割與打包
用對(duì)方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
質(zhì)量投訴超六成
汽車觀察(2016年3期)2016-02-28 13:16:26
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
主站蜘蛛池模板: 久久青青草原亚洲av无码| 国产精品一区在线观看你懂的| 久久香蕉欧美精品| 国产精品开放后亚洲| 国产成人精品在线1区| 综合人妻久久一区二区精品 | 91蜜芽尤物福利在线观看| 久久精品这里只有国产中文精品| 久久夜夜视频| 91偷拍一区| 手机永久AV在线播放| 综合成人国产| 欧美精品啪啪一区二区三区| 18禁不卡免费网站| 亚洲最黄视频| 中文字幕 91| 亚洲人成人无码www| 伊人成人在线视频| 成年网址网站在线观看| 日韩人妻无码制服丝袜视频| 71pao成人国产永久免费视频| 国产青青草视频| 国产簧片免费在线播放| 国产va视频| 国产美女丝袜高潮| 欧美黄网在线| a级免费视频| 亚洲熟女中文字幕男人总站| 最新日本中文字幕| 免费国产高清视频| 欧美日韩精品在线播放| 国产黄在线免费观看| 97国产在线视频| 97久久免费视频| 免费在线成人网| 伊人色婷婷| 亚洲国产欧美目韩成人综合| 亚洲成人在线免费| 找国产毛片看| 美美女高清毛片视频免费观看| 精品视频福利| 1024国产在线| 欧美不卡视频一区发布| 日本亚洲国产一区二区三区| 一级一级一片免费| 精品一区国产精品| 狠狠做深爱婷婷久久一区| 亚洲性影院| 国产精品免费露脸视频| 亚洲全网成人资源在线观看| 亚洲第一成年人网站| 国产好痛疼轻点好爽的视频| 国产精品久久精品| 国产极品美女在线播放| 精久久久久无码区中文字幕| 久精品色妇丰满人妻| 第一页亚洲| 色国产视频| 婷婷综合缴情亚洲五月伊| 精品无码国产一区二区三区AV| 亚洲精品无码专区在线观看| 欧美另类一区| 在线观看国产网址你懂的| 亚洲人成网站色7777| 亚洲成人网在线播放| 91美女在线| 尤物特级无码毛片免费| 精品国产自| 欧美97色| 国产精品亚洲专区一区| 中文精品久久久久国产网址 | 在线播放真实国产乱子伦| 午夜日韩久久影院| 爆操波多野结衣| 亚洲国产欧美自拍| 在线无码九区| 在线国产欧美| 午夜电影在线观看国产1区| 丁香六月综合网| 色网站免费在线观看| 丝袜国产一区| 欧美成人h精品网站|