王宏偉, 陳瑜瀟

(1.新疆大學電氣工程學院, 烏魯木齊 830047; 2.大連理工大學控制科學與工程學院， 大連 116024)

在工業(yè)實際生產過程中，受傳感器采樣頻率、物理設備、人工采樣、機械振動、網(wǎng)絡傳輸?shù)纫蛩氐挠绊?，使得系統(tǒng)采樣呈現(xiàn)輸入信號的刷新時間和輸出信號的采樣時間是不相同的, 此類系統(tǒng)是多采樣率數(shù)據(jù)系統(tǒng)。對于多采樣率線性數(shù)據(jù)的線性系統(tǒng)，利用變量提升技術[1-5]、多項式轉換技術[6-10]。結合多信息最小二乘類辨識算法，在不同模型結構下，如ARMAX(auto-regressive moving average exogenous)、Box-Jenkins、輸出誤差類模型等，研究者們創(chuàng)新性地提出了很多辨識算法，取得了很好的辨識建模效果。但是，實際工業(yè)系統(tǒng)往往呈現(xiàn)非線性，很難使用線性系統(tǒng)的辨識方法來辨識多采樣率的非線性系統(tǒng)。本文以雙采樣率數(shù)據(jù)的非線性系統(tǒng)為辨識對象開展研究，圖1給出此類系統(tǒng)的結構圖。

對于雙率采樣的Hammerstein非線性系統(tǒng)，文獻[11]提出了變量提升技術和最小二乘算法辨識方法，但是利用變量提升技術后的系統(tǒng)往往向量維數(shù)較大，對算法的收斂性沒有分析；文獻[12]針對含有有色噪聲的多采樣率數(shù)據(jù)非線性系統(tǒng)，提出基于遞階原理的隨機梯度辨識算法解決含有有色噪聲的多采樣率非線性系統(tǒng)辨識問題，算法收斂性也進行了探討；文獻[13]采用多新息隨機梯度辨識方法分別對非均勻采樣數(shù)據(jù)Hammerstein-Wiener系統(tǒng)和Wiener系統(tǒng)進行了建模研究；針對采用提升技術往往使系統(tǒng)變量、參數(shù)個數(shù)增加，模型結構復雜的問題，文獻[14]基于遞階原理將含有提升變量的系統(tǒng)模型辨識轉化為若干個子系統(tǒng)的辨識，這樣算法簡單、計算量小，效率高；針對含飽和死區(qū)特性的Hammerstein系統(tǒng)，文獻[15]提出了借助遞推思想和輔助模型思想估計內部未知變量，在此基礎上利用梯度辨識法進行參數(shù)進行估計。文獻[16-17]通過模糊建模方法解決了非均勻采樣的復雜非線性系統(tǒng)辨識問題。采用模糊聚類解決模糊模型前提結構和參數(shù)問題，利用遞階最小二乘算法解決模糊模型后件參數(shù)估計，同時對算法收斂性進行了研究。文獻[18]通過新型粒子濾波器而非輔助模型來估計不可測輸出量，利用梯度迭代算法進行系統(tǒng)辨識，仿真效果表明該方法與經典輔助模型算法相比更有效。文獻[19]將噪聲分布利用t分布表示，將未知時間延遲作為潛在變量，利用條件期望最大化對時間延遲未知且具有異常值的雙率輸入非線性方程誤差系統(tǒng)進行了參數(shù)辨識，仿真效果表明該方法辨識結果更精確。文獻[20]令同一被估計參數(shù)利用前三次估計值迭代得到新估計值，與傳統(tǒng)基于輔助模型的隨機梯度算法相比提高了收斂速率，且采用改進的卡爾曼濾波方法提高了辨識精度。

針對雙率采樣的Hammerstein非線性系統(tǒng)的辨識，一般都是針對模塊化的Hammerstein模型、Wiener模型、Hammerstein-Wiener模型、Wiener-Hammerstein模型展開的，其非線性模型一般是由靜態(tài)非線性環(huán)節(jié)和線性動態(tài)環(huán)節(jié)串聯(lián)組成的，其中非線性環(huán)節(jié)一般是由非線性多項式組成的，但實際工業(yè)系統(tǒng)中也有其他靜態(tài)非線性特性，例如滯后特性，死區(qū)特性、間隙特性、繼電器特性、飽和特性等，其中飽和特性是實際控制中經常遇到的，例如化工過程的電磁閥門輸入輸出特性、控制器放大器的飽和特性等。這些非線性系統(tǒng)的辨識遇到的困難主要為：①被辨識參數(shù)間存在交叉耦合和交叉乘積項；②含有很多未知的中間變量；③含有采集不到的缺失數(shù)據(jù)。本文針對含有飽和特性的非線性系統(tǒng)的辨識問題，基于雙采樣率數(shù)據(jù)開展了辨識建模的研究。

1 含有飽和特性的雙采樣率的Hammerstein系統(tǒng)模型

圖1 雙采樣率數(shù)據(jù)的Hammerstein系統(tǒng)結構Fig.1 Structure of Hammerstein system with dual-rate sampling data

對于圖1的Hammerstein非線性系統(tǒng)，其中線性動態(tài)環(huán)節(jié)為

y(t)=x(t)+v(t)

(1)

(2)

式(1)中：A(z)和B(z)是單位后移算子z-1[z-1y(t)=y(t-1)]的多項式，它們分別為A(z)=1+a1z-1+a2z-2+…+anaz-na、B(z)=b1z-1+b2z-2+…+bnbz-nb，其中ai和bi是多項式中的參數(shù)。本文中取b1=1。

圖1中，靜態(tài)非線性環(huán)節(jié)f(·)是飽和特性非線性函數(shù)，其輸入輸出關系如圖2所示。

圖2 飽和特性非線性函數(shù)Fig.2 Function of saturation nonlinearity

飽和特性非線性函數(shù)可以表示為

(3)

式(3)中：k1、k2分別為正、負線性部分的斜率；k1d1、k2d2分別為正、負線性部分的上、下限值。式(3)是分段線性函數(shù)，為了計算方便，可以將式(3)寫為

k1u(t)h[u(t)-d1]-k1d1h[u(t)]h×

[d1-u(t)]+k2u(t)h[-u(t)]+

k2d2h[-u(t)]-k2u(t)h[d2-u(t)]-

k2d2h[-u(t)]h[u(t)-d2]

(4)

式(4)中：h(t)是開關函數(shù)，即

(5)

對于式(5)，可以轉變?yōu)橐恍┓蔷€性基函數(shù)fj[u(t)]的展開，即

c4f4[u(t)]+c5f5[u(t)]+c6f6[u(t)]+

(6)

式(6)中：c1=k1；c2=k1d1；c3=k1；c4=k2；c5=k2d2；c6=k2；c7=k2；c8=k2d2；基函數(shù)為f1[u(t)]=u(t)h[u(t)]，f2[u(t)]=h[u(t)]，f3[u(t)]=-u(t)h[u(t)-d1]，f4[u(t)]=-h[u(t)]h[d1-u(t)]，f5[u(t)]=u(t)h[-u(t)]，f6[u(t)]=h[-u(t)]，f7[u(t)]=-u(t)h[d2-u(t)]，f8[u(t)]=-h[-u(t)]h[u(t)-d2]。

當然，式(6)又可以寫成向量表達式：

(7)

式(7)中：f[u(t)]={f1[u(t)],f2[u(t)],f3[u(t)],f4[u(t)],f5[u(t)],f6[u(t)],f7[u(t)],f8[u(t)]}T∈R8,c=[c1,c2,c3,c4,c5,c6,c7,c8]T∈R8。

定義輸入信息向量φ(t)和參數(shù)向量θ分別為

φ(t)={-x(t-1),-x(t-2),…,-x(t-na)

這樣根據(jù)式(1)、式(2)就有：

y(t)=φT(t)θ+v(t)

(8)

辨識的目的：利用雙率采樣數(shù)據(jù){u(k),y(kq)}，k=1,2,3，…,對靜態(tài)非線性模塊f(·)和動態(tài)線性模塊G(z)中的參數(shù)進行辨識。為此，用kq代替t就有：

y(kq)=φT(kq)θ+v(kq)

(9)

(2)非線性基函數(shù)fj[u(kq+i)],j=1,2,…,8,i=0,1,…,q-1也是未知的。為此，借助于輔助模型的輸出作為非線性基函數(shù)的估計。

2 基于輔助模型原理的遞推辨識

為了實現(xiàn)參數(shù)辨識，定義準則函數(shù)為

(10)

(11)

(12)

(13)

定義變量：

(14)

(15)

i=0,1,…,q-1

(16)

i=0,1,…,q-1

(17)

i=0,1,…,q-1

(18)

i=0,1,…,q-1

(19)

(20)

u(kq+i)],i=0,1,…,q-1

(21)

i=0,1,…,q-1

(22)

i=0,1,…,q-1

(23)

(24)

(25)

i=0,1,…,q-1

(26)

定義多個輔助模型后，基于輔助模型的遞推最小二乘辨識算法如下：

(27)

P(kq)=P(kq-q)-

(28)

i=0,1,…,q-1

(29)

i=0,1,…,q-1

(30)

i=0,1,…,q-1

(31)

i=0,1,…,q-1

(32)

下面將整個辨識算法完整總結如下：

(33)

(34)

P(kq)=P(kq-q)-

(35)

(36)

i=0,1,…,q-1

(37)

i=0,1,…,q-1

(38)

i=0,1,…,q-1

(39)

i=0,1,…,q-1

(40)

(41)

(42)

(43)

u(kq+i)]

(44)

(45)

(46)

u(kq+i)]

(47)

(48)

(49)

i=0,1,…,q-1

(50)

i=0,1,…,q-1

(51)

下面將基于雙率采樣數(shù)據(jù)的Hammerstein系統(tǒng)辨識過程總結如下：

(3)根據(jù)式(35)計算P(kq)。

(6)令k=k+1, 如果k≤N轉到步驟(2)。

圖3給出上述參數(shù)辨識算法的流程框圖。

圖3 參數(shù)辨識算法的流程框圖Fig.3 Flow chart of parameter identification algorithm

3 仿真算例

考慮下列含飽和特性的Hammerstein系統(tǒng)：

(52)

飽和特性非線性函數(shù)為

(53)

式中：A(z)=1+0.4z-1+0.3z-2，B(z)=1+0.5z-1+0.6z-2；v(k)分別是零均值，方差為σ2=0.0152、0.0252、0.0452的白噪聲；飽和特性非線性函數(shù)的參數(shù)為：k1=0.5,k2=0.3,d1=0.5,d2=-0.8。取雙率系數(shù)q=2，p0=1×106。

對于仿真實例的辨識模型，需要辨識的參數(shù)向量θ=[a1a2b1b2c1c2c3c4c5c6c7c8]T=[0.4 0.3 0.5 0.6 0.5 0.25 0.5 0.25 0.3 -0.24 0.3 -0.24]T。

圖4～圖11、表1～表6給出了辨識的結果，可以得到如下結論。

(1)在不同噪信比下參數(shù)估計誤差都隨著t的增加而減小，而且估計誤差減小速率很快，辨識精度較高。

(2)在不同噪信比下，參數(shù)估計相對誤差分別為1.148 5%、1.244 2%、3.083 1%，誤差變化不大，說明在不同噪信比下，在雙率采樣數(shù)據(jù)下，提出的辨識算法都能夠滿足參數(shù)估計的要求。

(3)在不同噪信比下，噪信比越大，參數(shù)估計精度越低，且估計值逼近真值的速率越低。

圖4 相對參數(shù)估計誤差曲線Fig.4 Relative parameter estimation error curve

4 結論

基于輔助模型思想, 提出輔助模型的最小二乘算法辨識雙采樣率的含飽和特性的Hammerstein非線性系統(tǒng)，仿真結果驗證了該算法的有效性。Hammerstein非線性系統(tǒng)雖然已經有了大量的研究基礎，但是其在雙率采樣條件下還有非常廣泛的、深入的研究空間，例如滯后非線性特性，死區(qū)非線性特性，繼電器非線性特性等，特別是繼電器非線性特性的研究較少，因此這方面后續(xù)的研究是非常有意義。

圖5 參數(shù)a、b估計值曲線Fig.5 Parameter estimation curve of a,b

圖6 σ2=0.0152時參數(shù)c估計值曲線Fig.6 Parameter estimation curve of c while σ2=0.0152

圖7 圖6局部放大圖 Fig.7 Detailed view from fig.6

圖8 σ2=0.0252時參數(shù)c估計值曲線Fig.8 Parameter estimation curve of c while σ2=0.0252

圖9 圖8局部放大Fig.9 Detailed view from fig.8

圖10 σ2=0.0452時參數(shù)c估計值曲線Fig.10 Parameter estimation curve of c while σ2=0.0452

圖11 圖10局部放大Fig.11 Detailed view from fig.10

表1 辨識得到的參數(shù)a、b及相對估計誤差(σ2=0.0152，δns=6.705 9%)Table 1 Parameter estimation of a,b and relative estimation error(σ2=0.0152，δns=6.705 9%)

表2 辨識得到的參數(shù)c(σ2=0.0152，δns=6.705 9%)Table 2 Parameter estimation of c (σ2=0.0152，δns=6.705 9%)

表3 辨識得到的參數(shù)a、b相對及估計誤差(σ2=0.0252，δns=11.216 6%)Table 3 Parameter estimation of a,band relative estimation error(σ2=0.0252，δns=11.216 6%)

表4 辨識得到的參數(shù)c(σ2=0.0252，δns=11.216 6%)Table 4 Parameter estimation of c (σ2=0.0252，δns=11.216 6%)

表6 辨識得到的參數(shù)c(σ2=0.0452，δns=20.117 7%)Table 6 Parameter estimation of c (σ2=0.0452，δns=20.117 7%)