熱電偶動(dòng)態(tài)響應(yīng)的帶外部輸入自回歸模型*

金敏俊,李文軍,鄭永軍,曾九孫

(中國計(jì)量大學(xué)計(jì)量測試工程學(xué)院,杭州 310018)

工業(yè)領(lǐng)域中熱電偶被廣泛應(yīng)用在各種溫度測量系統(tǒng)中[1]。隨著測量技術(shù)的發(fā)展,熱電偶也越來越多地應(yīng)用在各種介質(zhì)的瞬態(tài)溫度測量中[2-3],尤其是流體介質(zhì)瞬態(tài)溫度的測量[4]。當(dāng)熱電偶被用來測量瞬態(tài)溫度時(shí),需要對(duì)熱電偶進(jìn)行動(dòng)態(tài)響應(yīng)特性評(píng)價(jià)和校準(zhǔn)。隨著動(dòng)態(tài)測量技術(shù)的進(jìn)步,產(chǎn)生了多種評(píng)價(jià)和校準(zhǔn)方法。

近年CANUTO E等分別從理論和實(shí)驗(yàn)兩個(gè)角度出發(fā)[5],提出了動(dòng)態(tài)相對(duì)校準(zhǔn)的方法,討論了多個(gè)溫度傳感器組合測量和校準(zhǔn),包括產(chǎn)生精確溫度梯度問題、保持均勻穩(wěn)定溫度場問題。趙學(xué)敏等采用火焰為溫度激勵(lì)源[6],對(duì)氣體燃燒環(huán)境下的熱電偶動(dòng)態(tài)響應(yīng)特性進(jìn)行了分析。DINIZ A.C.G.C.A.等提出一種轉(zhuǎn)動(dòng)式校準(zhǔn)裝置和校準(zhǔn)方法[7],提供氣流環(huán)境下的溫度階躍,對(duì)熱電偶進(jìn)行動(dòng)態(tài)校準(zhǔn)并做不確定度分析,該方法計(jì)算了溫度對(duì)時(shí)間二階導(dǎo)數(shù)的極值,以極值對(duì)應(yīng)的時(shí)間點(diǎn)為階躍發(fā)生的時(shí)間。

在采用上述方法考察熱電偶動(dòng)態(tài)響應(yīng)時(shí),都需要確定溫度激勵(lì)產(chǎn)生的時(shí)間起點(diǎn)。用輔助測量手段確定時(shí)間起點(diǎn),會(huì)引入熱電偶輸入端的誤差。此外,目前的激勵(lì)方法中,施加的激勵(lì)溫度被視為一個(gè)穩(wěn)定的準(zhǔn)確值,但是如果采用流體比如氣體作為激勵(lì)介質(zhì),氣流溫度在中心值附近存在波動(dòng)。因此,當(dāng)以溫度激勵(lì)作為熱電偶的輸入,以熱電偶響應(yīng)為輸出,對(duì)熱電偶動(dòng)態(tài)響應(yīng)進(jìn)行統(tǒng)計(jì)建模時(shí)[8],由于輸入帶有誤差[9],所構(gòu)成的系統(tǒng)是一種變量帶誤差EIV(Error in Variables)系統(tǒng)[10],這種系統(tǒng)的辨識(shí)問題屬于變量帶誤差辨識(shí)問題[11]。本文從統(tǒng)計(jì)建模角度,建立了雙氣流環(huán)境下的熱電偶動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)系統(tǒng),考察熱電偶在高低溫兩種氣流之間擺動(dòng)時(shí)的動(dòng)態(tài)響應(yīng)。給出了描述熱電偶動(dòng)態(tài)響應(yīng)的狀態(tài)空間方程,采用一種帶外部輸入自回歸模型ARX(Autoregressive with exogenous,)對(duì)熱電偶動(dòng)態(tài)響應(yīng)過程進(jìn)行辨識(shí)[12]。以一種露端式鎳鉻鎳硅(Nickel-chromium/Nickel-silicon)熱電偶為實(shí)驗(yàn)對(duì)象,給出了一個(gè)算例。

1 熱電偶動(dòng)態(tài)響應(yīng)的實(shí)驗(yàn)方法及實(shí)驗(yàn)系統(tǒng)

1.1 熱電偶動(dòng)態(tài)響應(yīng)的實(shí)驗(yàn)方法及影響因素

熱電偶動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)方法主要包括熱風(fēng)洞法、電加熱法、激波管法、投入法和激光加熱法等。這些方法都是對(duì)熱電偶產(chǎn)生一個(gè)近似于理想階躍的溫度激勵(lì),通過測量熱電偶溫度響應(yīng)隨時(shí)間的變化,然后評(píng)價(jià)熱電偶動(dòng)態(tài)響應(yīng)能力。由于熱電偶測量接點(diǎn)與介質(zhì)之間的換熱受到多種換熱因素的影響,熱電偶動(dòng)態(tài)響應(yīng)也相應(yīng)地受到這些換熱因素的影響[13]。

以熱風(fēng)洞法為例分析,熱風(fēng)洞法產(chǎn)生溫度激勵(lì)的介質(zhì)為氣體,影響熱電偶與氣體介質(zhì)之間換熱的主要因素包括:對(duì)流以及流動(dòng)狀態(tài)、介質(zhì)物理性質(zhì)以及換熱面形狀、接觸面幾何尺寸和相對(duì)位置[14]。如果溫度激勵(lì)幅值較大,則熱電偶熱端指向冷端的熱傳導(dǎo)也不可忽略。同時(shí),熱電偶與氣體之間輻射換熱、與管道內(nèi)壁面之間也不可忽略。

如上所述,在氣流環(huán)境下存在較多因素影響熱電偶的動(dòng)態(tài)響應(yīng)。如果試圖從物理過程建模,對(duì)這些因素進(jìn)行選取和舍棄都比較困難。于是從實(shí)驗(yàn)數(shù)據(jù)角度建立模型成為一種合適的選擇,即利用統(tǒng)計(jì)建模和數(shù)據(jù)分析來估計(jì)影響因素、以及影響因素的顯著性?？紤]到氣體是測量動(dòng)態(tài)溫度的常見被測介質(zhì),因此采用常見的空氣為激勵(lì)介質(zhì),建立了一種空氣環(huán)境下的熱電偶動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)系統(tǒng)。

1.2 熱電偶動(dòng)態(tài)響應(yīng)的實(shí)驗(yàn)系統(tǒng)

根據(jù)上述分析,建立了一種氣流環(huán)境下的熱電偶動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)系統(tǒng)。為了避免管道壁面輻射引起的影響,以熱風(fēng)槍氣流出口為溫度激勵(lì)點(diǎn)。熱風(fēng)槍出口段加裝了整流段,以使氣流流動(dòng)均勻。熱電偶定位點(diǎn)如圖1所示。

圖1 熱電偶在管道出口的換熱模型

實(shí)驗(yàn)系統(tǒng)主要包括計(jì)算機(jī)、數(shù)據(jù)采集卡、可編程控制器、機(jī)械擺動(dòng)驅(qū)動(dòng)機(jī)構(gòu)和熱風(fēng)槍。其示意圖如圖2所示。

圖2 熱電偶動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)系統(tǒng)示意圖

實(shí)驗(yàn)原理是利用熱風(fēng)槍A、熱風(fēng)槍B分別產(chǎn)生高低溫氣流,用可編程控制器控制機(jī)械擺動(dòng)驅(qū)動(dòng)機(jī)構(gòu),控制并驅(qū)動(dòng)熱電偶測量接點(diǎn)在兩種氣流之間擺動(dòng),由數(shù)據(jù)采集卡采集熱電偶輸出,由可編程控制器返回?zé)犭娕歼M(jìn)入氣流的各時(shí)間點(diǎn)。

實(shí)驗(yàn)中,數(shù)據(jù)采集卡采用MCC USB1616HS高速數(shù)據(jù)采集卡,可編程控制器采用OMRON CPM1A型,熱風(fēng)槍采用QUICK2008型。

1.3 鎳鉻鎳硅熱電偶的動(dòng)態(tài)響應(yīng)數(shù)據(jù)集

以一種K型露端式鎳鉻鎳硅熱電偶為實(shí)驗(yàn)對(duì)象,利用圖2所示實(shí)驗(yàn)系統(tǒng),按照表1給出的實(shí)驗(yàn)配置參數(shù)進(jìn)行實(shí)驗(yàn)。

表1 動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)配置參數(shù)

按照上述配置參數(shù)進(jìn)行實(shí)驗(yàn),獲得了三支K型露端式鎳鉻鎳硅熱電偶的響應(yīng)數(shù)據(jù)集。圖3給出了低溫氣流為326 K且高溫氣流為515 K時(shí),這三支熱電偶的輸入和輸出。圖4給出了低溫氣流溫度為326 K且高溫氣流溫度為752 K時(shí),這三支熱電偶所對(duì)應(yīng)的輸入和輸出。

圖3 氣流溫度326 K和515 K時(shí)熱電偶的響應(yīng)

圖4 氣流溫度326 K和752 K時(shí)熱電偶的響應(yīng)

通過實(shí)驗(yàn)所獲得的響應(yīng)數(shù)據(jù)都帶有測量誤差,如果分別用來和ym(t)表示熱電偶的輸入和輸出,則熱電偶動(dòng)態(tài)響應(yīng)過程如圖5所示。

圖5 輸入和輸出帶有誤差的過程

圖5中,um(t)為熱電偶輸入的測量值并含有干擾eu(t),ym(t)為熱電偶輸出的測量值并含有干擾ey(t),yz(t)為作用在熱電偶輸出上的加性干擾。干擾ey(t)和yz(t)可以合并處理,但eu(t)的存在,使得上述熱電偶動(dòng)態(tài)響應(yīng)過程構(gòu)成一類變量帶誤差回歸問題[15]。因此,把熱電偶視為一種帶有隨機(jī)誤差的確定性系統(tǒng),建立其狀態(tài)空間模型,利用帶外部輸入自回歸方法進(jìn)行辨識(shí)。為了表達(dá)簡便,本文以下的描述中熱電偶的輸入um(t)和輸出ym(t)仍分別用u(t)和表示。

2 熱電偶動(dòng)態(tài)響應(yīng)的自回歸模型

2.1 熱電偶的連續(xù)時(shí)間狀態(tài)空間模型

熱電偶動(dòng)態(tài)響應(yīng)的實(shí)驗(yàn),是由激勵(lì)介質(zhì)給予熱電偶溫度激勵(lì),這種激勵(lì)包括正向溫度階躍和負(fù)向溫度階躍。熱電偶受到溫度激勵(lì)后,響應(yīng)隨著時(shí)間連續(xù)變化。把熱電偶視為一個(gè)單輸入單輸出連續(xù)時(shí)間系統(tǒng),其狀態(tài)空間方程可以描述為:

(1)

y(t)=cx(t)+du(t)

(2)

式(1)為熱電偶狀態(tài)方程,式(2)為熱電偶輸出方程。式中,x(t)為熱電偶的狀態(tài)變量,y(t)為熱電偶輸出,u(t)為熱電偶輸入,A∈Rn×n為系數(shù)矩陣,b∈R1×n、c∈R1×n、d∈R1×n都是系數(shù)矩陣。

把式(1)和式(2)寫成解的形式,其描述為[25]:

(3)

(4)

在熱電偶動(dòng)態(tài)響應(yīng)的實(shí)驗(yàn)中,需要通過采樣方法獲得熱電偶的輸入和輸出?？紤]到施加的溫度激勵(lì)通常為正向階躍或者負(fù)向階躍,即輸入近似為方波信號(hào),采用零階保持器,利用階躍響應(yīng)不變變換把連續(xù)系統(tǒng)轉(zhuǎn)化為離散系統(tǒng)。設(shè)采樣周期為T,離散時(shí)間變量為k,采用零階保持器,即:

u(t)=u(kT)

(5)

kT≤t≤(k+1)T

(6)

再定義:

x(k)=x(kT)

(7)

u(k)=u(kT)

(8)

y(k)=y(kT)

(9)

這時(shí),式(3)和式(4)所表示的熱電偶連續(xù)時(shí)間狀態(tài)空間模型轉(zhuǎn)化為離散時(shí)間狀態(tài)空間模型:

x(k+1)=Gx(k)+fu(k)

(10)

y(k)=cx(k)+du(k)

(11)

式中:G=eAt,f=A-1[G-I]b,I為單位陣。

式(10)和式(11)給出了在確定的采樣周期下,采用零階保持器時(shí),熱電偶連續(xù)時(shí)間狀態(tài)空間模型所對(duì)應(yīng)的離散時(shí)間狀態(tài)空間模型的形式。在熱電偶動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)中,需要通過采樣獲得離散的輸入輸出數(shù)據(jù),因此使用離散時(shí)間系統(tǒng)建立模型。

2.2 熱電偶的離散時(shí)間狀態(tài)空間模型

把熱電偶響應(yīng)過程視為離散時(shí)間系統(tǒng),其中離散時(shí)間變量用k表示,這時(shí)熱電偶的狀態(tài)空間模型可以表示為:

x(k+1)=Ax(k)+bu(k)

(12)

y(k)=cx(k)+du(k)

(13)

式(12)和式(13)中,x(k)為熱電偶的狀態(tài)變量,y(k)為熱電偶的輸出序列,u(k)為熱電偶的輸入序列,A∈Rn×n為常系數(shù)矩陣,b∈R1×n、c∈R1×n、d∈R1×n都是系數(shù)矩陣。

引入移位算子,即單位前移算子z:

zy(k)=y(k+1)

(14)

以及單位后移算子z-1:

z-1y(k)=y(k-1)

(15)

式(12)和式(13)可以表示為:

zx(k)=Ax(k)+bu(k)

(16)

y(k)=cx(k)+du(k)

(17)

把式(16)、式(17)表示為解的形式,有:

x(k)=(zI-A)-1bu(k)

(18)

y(k)=[c(zI-A)-1b+ud(k)

(19)

式(18)和式(19)是描述熱電偶動(dòng)態(tài)響應(yīng)過程的離散時(shí)間系統(tǒng)。由于輸入和輸出都存在測量誤差,熱電偶動(dòng)態(tài)響應(yīng)過程可用隨機(jī)擾動(dòng)與確定性模型結(jié)合的模型來描述,即利用隨機(jī)系統(tǒng)模型對(duì)離散時(shí)間隨機(jī)信號(hào)的統(tǒng)計(jì)性質(zhì)進(jìn)行評(píng)估。

2.3 熱電偶的變量帶誤差自回歸模型

對(duì)于帶隨機(jī)信號(hào)的確定性模型,其典型結(jié)構(gòu)有自回歸模型AR(Autoregressive model)、帶外部輸入自回歸模型ARX(Autoregressive model with exogenous input)、帶外部輸入自回歸滑動(dòng)平均ARMAX(Autoregressive Moving Average model with exogenous input)、輸出誤差模型OE(Out Error model)和博克斯-詹金斯模型BJ(Box-Jenkins model)等多種。對(duì)于離散時(shí)間系統(tǒng),用多項(xiàng)式表示這幾種模型的一般表達(dá)式為:

(20)

式(20)中,A(z),B(z),C(z),D(z)和F(z)是用移位算子z表示的多項(xiàng)式,ui表示第i個(gè)輸入,nu表示輸入數(shù)目,nki表示第i個(gè)輸入的時(shí)延,n(k)為隨機(jī)干擾。

在上述幾種模型結(jié)構(gòu)中,ARX模型在工業(yè)中得到廣泛應(yīng)用,可用于變量帶誤差問題的處理[16-17]。ARX模型的一般表達(dá)式為:

(21)

對(duì)于熱電偶動(dòng)態(tài)響應(yīng)過程,如果忽略熱電偶進(jìn)出氣流的移動(dòng)速度等次要影響因素,式(21)可簡化為單輸入單輸出ARX模型[18]:

A(z)y(k)=B(z)u(k-nk)+n(k)

(22)

即:

(23)

式(23)中,y(k)為熱電偶輸出序列,A(z)和B(z)是用移位算子z表示的多項(xiàng)式,u(k)為熱電偶輸入序列,nk表示輸入時(shí)延,n(k)表示干擾序列。定義na為輸出時(shí)延長度,nb為輸入時(shí)延長度,則多項(xiàng)式A(z)可表示為常系數(shù)多項(xiàng)式:

A(z)=1+a1z-1+a2z-2+…+anaz-na

(24)

多項(xiàng)式B(z)可表示為:

B(z)=b1z-1+b2z-2+…+bnbz-nb

(25)

對(duì)應(yīng)的模型示意圖如圖6。

圖6 帶外部輸入自回歸模型

圖7 變量帶誤差系統(tǒng)的帶外部輸入自回歸模型

考慮熱電偶輸入u(k)帶有干擾,即變量帶誤差情形時(shí),式(23)應(yīng)表示為:

(26)

u(k)=u0(k)+eu(k)

(27)

對(duì)應(yīng)的ARX模型如圖7所示。

對(duì)于圖7所示的模型,當(dāng)用線性模型時(shí),定義參數(shù)向量為:

θ=[a1,a2,…,ana,b1,b2,…,bnb]T

(28)

定義數(shù)據(jù)向量為:

φ=[-y(k-1),-y(k-2),…,-y(k-na),u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)]T

(29)

并定義:

(30)

則由式(26)和式(27)可得:

(31)

(32)

式(31)包含熱電偶當(dāng)前輸出溫度y(k)值,以及輸出時(shí)延值、輸入當(dāng)前值和輸入時(shí)延值。把當(dāng)前輸出溫度y(k)的預(yù)報(bào)值記為yp(k),則由式(31)可構(gòu)造當(dāng)前值的一種預(yù)報(bào)值,其展開式可以表示為:

yp(k)=[-a1,-a2,…,-ana,b1,b2,…,bnb]× [y(k-1),y(k-2),…,y(k-na),u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)]T

(33)

式(33)中,時(shí)延值y(k-1),y(k-2),…,y(k-na)和u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)為回歸量,-a1,-a2,…,-ana,b1,b2,…,bnb為回歸量的系數(shù),表示回歸量的權(quán)重。

當(dāng)用非線性函數(shù)f預(yù)報(bào)y(k),記:

yp(k)=f[y(k-1),y(k-2),…,y(k-na),u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)]

(34)

式(34)構(gòu)成了預(yù)報(bào)yp(k)的非線性ARX模型。

式(33)和式(34)可以統(tǒng)一表示為圖8所示的模型結(jié)構(gòu),其中左側(cè)矩形框表示計(jì)算估計(jì)量所采用的回歸量,右側(cè)矩形框表示計(jì)算估計(jì)量所采用的非線性函數(shù)和線性函數(shù)。

圖8 ARX模型的回歸量和估計(jì)量

圖8也描述了預(yù)報(bào)yp(k)的步驟。其中第一個(gè)步驟是計(jì)算回歸量。對(duì)于不同類型的熱電偶,可以采用不同類型的回歸量。如果采用形如u(k)、u(k-1)、y(k-1)的形式,是線性形式。如果采用形如u(k-1)×y(k-1)的形式,是u(k-1)和y(k-1)所構(gòu)成的非線性形式。

第二個(gè)步驟是用回歸量預(yù)報(bào)估計(jì)量。對(duì)于非線性函數(shù),仍用φ表示第一個(gè)步驟中的回歸量,圖8中計(jì)算估計(jì)量的非線性函數(shù)的一般表達(dá)式表示如下[19]:

(35)

θ=[αk,βk,γk]k=1,2,…,d

(36)

式(35)和式(36)中,θ為參數(shù)向量,φ為向量形式的回歸量,κ是非線性函數(shù),αk、βk和γk為參數(shù)。本文算例中非線性函數(shù)κ采用了工業(yè)上常用的 Sigmoid 函數(shù),函數(shù)形式為[20]:

κ(x)=1/(1+e-x)

(37)

2.4 準(zhǔn)則函數(shù)和評(píng)價(jià)指標(biāo)

采用式(31)所表示的線性模型和式(35)所表示的非線性模型,能夠?qū)犭娕紕?dòng)態(tài)響應(yīng)數(shù)據(jù)做有限數(shù)據(jù)長度的離線辨識(shí)。對(duì)于線性模型或非線性模型,或者兩者的混合模型,y的單步預(yù)報(bào)值都可表示為[21]:

y(k)=f(θ,φ)+e(k)

(38)

式(38)中,θ為參數(shù)向量,φ為數(shù)據(jù)向量,且有:

φ=[y(k-1),y(k-2),…,y(k-na),u(k-nk),u(k-nk-1),…,u(k-nk-nb+1)]T

(39)

由式(39),單步預(yù)報(bào)誤差為:

e(k)=y(k)-f(θ,φ)

(40)

定義準(zhǔn)則函數(shù)為:

J(θ)=y(k)-f(θ,φ)2

(41)

再對(duì)準(zhǔn)則函數(shù)J(θ)極小化,有:

(42)

根據(jù)上述準(zhǔn)則函數(shù),利用最小二乘法計(jì)算參數(shù)θ。在收斂情形下,可獲得最優(yōu)的θ和對(duì)應(yīng)的模型。

在獲得模型后,用擬合度FIT對(duì)模型進(jìn)行評(píng)價(jià),具體形式如下:

(43)

同時(shí)計(jì)算模型的均方根誤差RMSE,具體形式如下:

(44)

3 仿真

3.1 仿真方案

為了驗(yàn)證上述模型和算法,先通過仿真的方法對(duì)一種已知熱電偶模型進(jìn)行辨識(shí)。仿真分為兩種辨識(shí)方法,第一種方法是用傳統(tǒng)的自回歸模型AR進(jìn)行辨識(shí)[22],第二種方法是用帶外部輸入自回歸模型ARX進(jìn)行辨識(shí)。

假設(shè)已知熱電偶模型用以下傳遞函數(shù)描述:

(45)

根據(jù)式(45),用偽二值序列作為輸入生成無噪聲輸出,并用數(shù)據(jù)集zr表示。數(shù)據(jù)集zr的時(shí)間步長為0.2,長度為501,如圖9所示。

圖9 已知熱電偶模型的無噪聲輸出和輸入

上述偽二值序列增加服從高斯分布N(0,0.012)的隨機(jī)噪聲,作為含噪聲的輸入,并由式(45)產(chǎn)生對(duì)應(yīng)的輸出,作為含噪聲的輸出,并用數(shù)據(jù)集zm表示,如圖10所示。

圖10 已知熱電偶模型的含噪聲輸出和輸入

對(duì)于圖9中的無噪聲輸出數(shù)據(jù)集zr,劃分為兩個(gè)部分,前334個(gè)數(shù)據(jù)為zr1,后167個(gè)數(shù)據(jù)為zr2;對(duì)于圖10中的含噪聲輸出數(shù)據(jù)集zm,同樣劃分為兩個(gè)部分,前334個(gè)數(shù)據(jù)為zm1,后167個(gè)數(shù)據(jù)為zm2。

3.2 自回歸模型

首先利用自回歸模型AR對(duì)含噪聲數(shù)據(jù)zm1做辨識(shí)。通過計(jì)算,確定AR模型的階數(shù)為4,對(duì)應(yīng)的回歸量為:y(k-1)、y(k-2)、y(k-3)、y(k-4)。圖11 和圖12給出了對(duì)應(yīng)的計(jì)算結(jié)果,其中圖11為自回歸模型AR輸出與數(shù)據(jù)zr1的比較,圖12為自回歸模型AR輸出與數(shù)據(jù)zr2的比較。

圖11 自回歸模型AR輸出與數(shù)據(jù)zr1

圖12 自回歸模型AR輸出與數(shù)據(jù)zr2

從圖11和圖12中可以看出,AR模型對(duì)zr1的擬合度FIT為53.22%,均方根誤差RMSE為0.075 0;AR模型對(duì)zr2的擬合度FIT為64.81%,均方根誤差RMSE為0.063 5。

3.3 帶外部輸入自回歸模型

在利用自回歸模型AR辨識(shí)的基礎(chǔ)上,再利用帶外部輸入自回歸模型ARX對(duì)含噪聲數(shù)據(jù)zm1做辨識(shí)。先通過線性ARX模型計(jì)算,得到na為5,nb為4,nk為1,確定回歸量為:y(k-1)、y(k-2)、y(k-3)、y(k-4)、y(k-5)以及u(k-1)、u(k-2)、u(k-3)、u(k-4)。再根據(jù)線性模型所確定的回歸量,采用線性和非線性ARX混合模型做辨識(shí),其中κ函數(shù)采用κ(x)=1/(1+e-x)函數(shù)。在選取估計(jì)量時(shí),線性部分和非線性部分使用相同的回歸量,混合后的形式為:

f(φ)=LT(φ-r)+d+κ(Q(φ-r))

(46)

式(46)中,φ為作為回歸量的數(shù)據(jù)向量,L為向量形式的線性回歸系數(shù),r為φ的均值,d為標(biāo)量形式的偏移量,Q為投影矩陣。上述ARX模型的計(jì)算,模型參數(shù)包含了L,r,d,Q以及表征κ函數(shù)的參數(shù)。圖13和圖14給出了對(duì)應(yīng)的計(jì)算結(jié)果。其中圖13 為ARX模型輸出與數(shù)據(jù)zr1的比較,圖14為ARX模型輸出與數(shù)據(jù)zr2的比較。

圖13 帶外部輸入自回歸模型ARX 輸出與數(shù)據(jù)zr1

圖14 帶外部輸入自回歸模型ARX 輸出與數(shù)據(jù)zr2

計(jì)算得到ARX模型對(duì)zr1的擬合度FIT為94.69%,均方根誤差RMSE為0.008 5,ARX模型對(duì)zr2的擬合度FIT為95.96%,均方根誤差RMSE為0.007 3。

比較AR模型和ARX模型的辨識(shí)結(jié)果可以看出,ARX模型的擬合度比AR模型高,并且ARX模型的均方根誤差RMSE比AR模型小。

4 算例

在仿真的基礎(chǔ)上,對(duì)圖3和圖4所示熱電偶動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)數(shù)據(jù)做ARX辨識(shí)。

把圖3中的實(shí)驗(yàn)數(shù)據(jù)做歸一化處理并平均,再把數(shù)據(jù)集劃分為兩個(gè)部分,取前3 500個(gè)數(shù)據(jù)z1做辨識(shí),取后1 750個(gè)數(shù)據(jù)z2做驗(yàn)證。對(duì)數(shù)據(jù)集z1,采用線性模型計(jì)算得到na為5,nb為5,nk為1,確定回歸量為:y(k-1)、y(k-2)、y(k-3)、y(k-4)、y(k-5)以及u(k-1)、u(k-2)、u(k-3)、u(k-4)、u(k-5)。根據(jù)線性模型所確定的回歸量,再采用線性和非線性ARX混合模型計(jì)算,其中κ函數(shù)采用κ(x)=1/(1+e-x)函數(shù)。圖15和圖16給出了對(duì)應(yīng)的計(jì)算結(jié)果,其中圖15為ARX模型輸出與數(shù)據(jù)z1的比較,圖16為ARX模型輸出與數(shù)據(jù)z2的比較。

圖15 帶外部輸入自回歸模型ARX輸出與數(shù)據(jù)z1

圖16 帶外部輸入自回歸模型ARX輸出與數(shù)據(jù)z2

計(jì)算得到ARX模型對(duì)z1的擬合度FIT為 89.07%,均方根誤差RMSE為0.032 2,ARX模型對(duì)z2的擬合度FIT為88.92%,均方根誤差RMSE為0.031 9。

類似地,對(duì)于圖4中的實(shí)驗(yàn)數(shù)據(jù)做歸一化處理并平均,把數(shù)據(jù)集劃分為兩個(gè)部分,取前3 500個(gè)數(shù)據(jù)z3做辨識(shí),取后1 750個(gè)數(shù)據(jù)z4做驗(yàn)證。對(duì)數(shù)據(jù)集z3,采用線性模型計(jì)算得到na為5,nb為4,nk為1,確定回歸量為:y(k-1)、y(k-2)、y(k-3)、y(k-4)、y(k-5)以及u(k-1)、u(k-2)、u(k-3)、u(k-4)。根據(jù)線性模型所確定的回歸量,再采用線性和非線性ARX混合模型計(jì)算,其中κ函數(shù)采用κ(x)=1/(1+e-x)函數(shù)。圖17和圖18給出了對(duì)應(yīng)的計(jì)算結(jié)果,其中圖17為ARX模型輸出與數(shù)據(jù)z3的比較,圖18為ARX模型輸出與數(shù)據(jù)z4的比較。

圖18 帶外部輸入自回歸模型ARX輸出與數(shù)據(jù)z4

圖17 帶外部輸入自回歸模型 ARX輸出與數(shù)據(jù)z3

計(jì)算得到ARX模型對(duì)z3的擬合度FIT為91.67%,均方根誤差RMSE為0.024 4,ARX模型對(duì)z4的擬合度FIT為91.64%,均方根誤差RMSE為0.025 0。

綜合圖15～圖18可以看出,對(duì)于圖3所示雙氣流之間溫度差較小的情形,以及圖4所示雙氣流之間溫度差較大的情形,ARX模型都得到了較好的擬合度FIT以及較小的均方根誤差RMSE。

5 結(jié)論

建立了一種雙氣流環(huán)境下熱電偶動(dòng)態(tài)響應(yīng)實(shí)驗(yàn)系統(tǒng),對(duì)熱電偶產(chǎn)生正負(fù)溫度激勵(lì)并采集了熱電偶響應(yīng)數(shù)據(jù)。建立了熱電偶動(dòng)態(tài)響應(yīng)的狀態(tài)空間方程,把熱電偶動(dòng)態(tài)響應(yīng)過程化簡為隨機(jī)加確定性系統(tǒng),采用帶外部輸入自回歸模型,利用線性函數(shù)和非線性函數(shù)逼近熱電偶動(dòng)態(tài)響應(yīng)函數(shù),處理了熱電偶動(dòng)態(tài)響應(yīng)回歸分析中存在的變量帶誤差問題。實(shí)驗(yàn)和計(jì)算結(jié)果表明,帶外部輸入自回歸模型適用于K型露端式熱電偶動(dòng)態(tài)響應(yīng)過程的回歸分析。上述工作為工業(yè)熱電偶動(dòng)態(tài)響應(yīng)能力評(píng)價(jià)提供了一種測試方法。這種方法還可應(yīng)用于熱電偶測量動(dòng)態(tài)溫度的遞推估計(jì)中。