方旭，劉金琨

(北京航空航天大學(xué) 自動(dòng)化科學(xué)與電氣工程學(xué)院，北京100191)

0 引言

目前，垂直起降飛行器和無人機(jī)(UAV)的研究主要集中在航跡規(guī)劃和跟蹤兩個(gè)方面，而對(duì)UAV的參數(shù)辨識(shí)研究并不多，特別是在飛行器參數(shù)集員辨識(shí)［1］方面缺少研究。在實(shí)際應(yīng)用中飛行器系統(tǒng)的噪聲統(tǒng)計(jì)特性很難確定，用統(tǒng)計(jì)類的辨識(shí)方法是行不通的，而集員辨識(shí)［2-4］卻適合這種情況。集員辨識(shí)是在未知但有界噪聲假設(shè)下進(jìn)行的辨識(shí)，完成的工作是在參數(shù)空間中找到一個(gè)與量測(cè)數(shù)據(jù)和已知噪聲界相容的可行解集，即參數(shù)空間的一個(gè)集合［5］，且集合中的每一個(gè)成員均為可行解。隨著樣本的增大，成員集范圍逐漸縮小［6-7］。差分進(jìn)化算法［8］在保留了遺傳算法優(yōu)點(diǎn)的同時(shí)，又避免了復(fù)雜的運(yùn)算。跟蹤微分器［9］通過積分提取微分信號(hào)，積分鏈?zhǔn)轿⒎制骺梢杂行У匾种圃肼?，而且可以直接估?jì)系統(tǒng)高階導(dǎo)數(shù)。

本文為了解決辨識(shí)模型回歸矩陣奇異條件下的參數(shù)集員辨識(shí)問題，基于差分進(jìn)化算法和非線性向量回歸［10］，建立了逼近參數(shù)向量與誤差向量的范數(shù)之間的函數(shù)關(guān)系模型，并用差分進(jìn)化算法優(yōu)化回歸模型參數(shù)，避免了模型參數(shù)選擇的盲目性。此方程的應(yīng)用不要求回歸矩陣非奇異，在奇異條件下也能進(jìn)行集員辨識(shí)，擴(kuò)大了適用范圍。首先利用跟蹤微分器由飛行器的位置和歐拉角信息估計(jì)其6個(gè)自由度的速度和加速度信息;然后，將整個(gè)系統(tǒng)分解為兩個(gè)相對(duì)獨(dú)立的子系統(tǒng)，簡(jiǎn)化辨識(shí)模型，利用差分進(jìn)化算法對(duì)飛行器參數(shù)可行集的中心進(jìn)行估計(jì)和建立回歸函數(shù);最后，根據(jù)得到的函數(shù)關(guān)系模型導(dǎo)出近似參數(shù)可行解集。

1 無人機(jī)系統(tǒng)動(dòng)態(tài)模型

圖1為四旋翼UAV受力圖，其中Oxyz為慣性坐標(biāo)系，飛行器有4個(gè)螺旋槳控制位置和歐拉角，l為飛行器半臂長(zhǎng)，F(xiàn)i(i=1，2，3，4)為螺旋槳推力形成的扭轉(zhuǎn)力矩。

圖1 四旋翼UAV結(jié)構(gòu)圖Fig.1 Structure of quadrotor UAV

利用歐拉-拉格朗日方法推導(dǎo)出UAV的動(dòng)力學(xué)方程，其簡(jiǎn)化表達(dá)式如下:

式中:x，y，z分別為位置信息;θ，φ，ψ 分別為俯仰角、滾轉(zhuǎn)角和偏航角;ai(i=1，…，10)為待辨識(shí)的參數(shù)，如表1所示，表中m為質(zhì)量，c為比例系數(shù)，κi(i=1，…，6)為阻力系數(shù)，Ii(i=1，2，3)為轉(zhuǎn)動(dòng)慣量;ui(i=1，2，3，4)為虛擬輸入，滿足:

表1 待辨識(shí)的系統(tǒng)參數(shù)Table 1 Unknown parameters for identification

2 跟蹤微分器和差分進(jìn)化算法

四旋翼UAV的位置和歐拉角可以通過傳感器來測(cè)量，設(shè)計(jì)了6個(gè)三階跟蹤微分器來分別獲得飛行器位置、歐拉角的速度和加速度數(shù)據(jù)，三階跟蹤微分器設(shè)計(jì)如下:

其中s3+c3s2+c2s+c1=0滿足Hurwitz穩(wěn)定判據(jù)，μ是極小的正數(shù)，輸入u是傳感器測(cè)量到的位置或歐拉角信息，y1，y2為輸出，方程(4)的Laplace變換為:

差分進(jìn)化算法是比較成熟的智能算法，本文不做詳細(xì)介紹，可參見文獻(xiàn)［11］。差分進(jìn)化算法的運(yùn)行參數(shù)主要有:變異因子F、交叉因子CR、群體規(guī)模Z和最大迭代次數(shù)T。其收斂速度主要由F和CR決定:F較小時(shí)，易陷入局部極值;F較大時(shí)，收斂速度會(huì)減慢。為此，可以采用下式線性調(diào)整［12］變異因子F。

式中:d為當(dāng)前進(jìn)化代數(shù);D為最大進(jìn)化代數(shù);Fmax和Fmin為選定的變異因子的最大、最小值。

CR越小，種群多樣性減小，過早收斂;CR越大，收斂速度越大，但收斂變慢。同樣，可以采用式(7)線性調(diào)整交叉因子CR:

為了保證算法的性能，CRmax和CRmin應(yīng)選取合理的值。隨著進(jìn)化代數(shù)的增加，F(xiàn)線性遞減，CR線性遞增，目的是希望改進(jìn)的差分進(jìn)化算法在搜索初期能夠保持種群的多樣性，到后期有較大的收斂速度。

3 飛行器系統(tǒng)分解和參數(shù)的中心估計(jì)

為了簡(jiǎn)化參數(shù)辨識(shí)的復(fù)雜度，將整個(gè)飛行器系統(tǒng)分解為兩個(gè)獨(dú)立的子系統(tǒng)，分別對(duì)其進(jìn)行參數(shù)辨識(shí)。這樣就成功地將一個(gè)高維的系統(tǒng)辨識(shí)問題降階分解到兩個(gè)低維的子系統(tǒng)上，避免了差分進(jìn)化算法在進(jìn)行高維優(yōu)化時(shí)容易陷入局部最優(yōu)的問題，同時(shí)也加快了系統(tǒng)的收斂速度。

兩個(gè)獨(dú)立的子系統(tǒng)(位置子系統(tǒng)(8)和姿態(tài)子系統(tǒng)(9))寫成含有噪聲的參數(shù)線性化回歸模型形式m=Nθ+e如下:

式中:e∈R3為未知但有界的噪聲;m∈R3為測(cè)量輸出。位置子系統(tǒng)(8)記為ma=Naθa+ea，其中Na∈R3×4，θa∈ R4;而姿態(tài)子系統(tǒng)(9)記為 mb=Nbθb+eb，其中Nb∈R3×6，θb∈R6。系統(tǒng)分解后就可以分別對(duì)子系統(tǒng)進(jìn)行集員辨識(shí)，下面以其中任意一個(gè)子系統(tǒng)m=Nθ+e為例介紹中心估計(jì)方法［13］。假設(shè)得到k組數(shù)據(jù)，記:

可得ms=Nsθ+es，集員辨識(shí)常用最小二乘法進(jìn)行參數(shù)的中心估計(jì)，但要求非奇異，限制了適用范圍。

因此，采用改進(jìn)差分進(jìn)化算法進(jìn)行中心估計(jì)，辨識(shí)誤差指標(biāo)選為:

圖2為集員辨識(shí)中心點(diǎn)估計(jì)流程圖。

圖2 集員辨識(shí)中心點(diǎn)估計(jì)流程圖Fig.2 Flow chart of set membership identification center estimation

4 基于差分算法和向量回歸的集員估計(jì)

以飛行器系統(tǒng)中任意一個(gè)子系統(tǒng)m=Nθ+e為例介紹集員辨識(shí)方法，與文獻(xiàn)［15-16］不同，本文采用差分進(jìn)化算法對(duì)回歸模型參數(shù)進(jìn)行優(yōu)化，并且給出適合四旋翼參數(shù)集員辨識(shí)的加權(quán)方法。e∈R3是有界的，滿足Ve∈R3是l2范數(shù)意義下有界的，即為對(duì)角加權(quán)矩陣，d∈R為誤差界限。由上可得與回歸模型、測(cè)量輸出、誤差界限相容的集合為:

式中:向量θ(i)(i=1，…，g)為參數(shù)空間中的點(diǎn)。根據(jù)樣本集S建立逼近誤差向量l2范數(shù)ε∈R與θ之間的函數(shù)關(guān)系的模型:

式中:η(·)為從低維參數(shù)空間θ∈Rp到高維特征空間Rq(q＞p)的非線性映射;w∈Rq為權(quán)重向量;b∈R為偏置量。為了得到w，b，依據(jù)非線性回歸結(jié)構(gòu)風(fēng)險(xiǎn)最小化原則，將其學(xué)習(xí)過程轉(zhuǎn)化為優(yōu)化問題如下:

式中:i=(1，…，l);Bi為松弛變量;γ為懲罰系數(shù)。利用拉格朗日方程求解最佳參數(shù):

ki∈R為拉格朗日乘子，求解最優(yōu)參數(shù)，得到下列矩陣方程:

式中:ξ =［ε1，…，εg］;lv=［1，…，1］;k=［k1，…，kg］;Ω∈Rg×g為矩陣中的元素，Ωij= η(θ(i))Tη(θ(j))(i，j=1，…，g)，根據(jù)梅塞條件存在核函數(shù)滿足式(14)。本文選取核 函數(shù)為σ2}(i，j=1，…，g)，其中 σ 為核函數(shù)寬度系數(shù)，解矩陣方程得到b，k，最終得到最優(yōu)方程解:

那么優(yōu)化問題式(12)變?yōu)?

同樣利用拉格朗日方程求解最佳參數(shù)，可得到最終模型形式:

因此集員估計(jì)可行解公式為:

假設(shè)得到Tk組采樣數(shù)據(jù)m，N，每組數(shù)據(jù)得到一個(gè)可行解:

如果真實(shí)參數(shù)的先驗(yàn)信息θ*∈Θ0，Θ0是參數(shù)空間的有界集合。那么，系統(tǒng)參數(shù)集員估計(jì)的可行解為Θ0與Hj的交集，即:

5 仿真實(shí)驗(yàn)與分析

改進(jìn)的差分進(jìn)化算法參數(shù)設(shè)置為:Fmax=0.9，F(xiàn)min=0.4，CRmax=0.9，CRmin=0.3，Z=80，D=500;跟蹤微分器參數(shù)設(shè)置為:c1=200，c2=140，c3=22，μ=0.01;飛行器輸入為:u1=10 sin(t)，u2=10 cos(t)，u3=10 sin(t)，u4=10 cos(t);系統(tǒng)物理參數(shù)為:a1=0.465 1，a2=0.051 2，a3=0.055 8，a4=0.051 2，a5=0.195 3，a6=0.033 2，a7=0.198 4，a8=0.031 7，a9=0.463 4，a10=0.013 1;最優(yōu)參數(shù)σ=0.025，γ=1e－4;根據(jù)測(cè)量輸出信號(hào)的變化范圍，取噪聲向量元素的界限為0.03。集員參數(shù)辨識(shí)結(jié)果如表2所示。

表2 集員參數(shù)辨識(shí)結(jié)果Table 2 Results for UAV set membership identification

圖3為系統(tǒng)6個(gè)自由度的幅值變化和系統(tǒng)噪聲。圖4為兩個(gè)子系統(tǒng)改進(jìn)差分進(jìn)化算法參數(shù)中心估計(jì)的辨識(shí)誤差指標(biāo)收斂過程。

圖3 6個(gè)自由度幅值變化和噪聲Fig.3 Amplitude change of 6 degress and noise

圖4 子系統(tǒng)辨識(shí)指標(biāo)Fig.4 Identification index of subsystem

從以上圖表可以發(fā)現(xiàn)，高階積分鏈?zhǔn)轿⒎制骺梢跃_地跟蹤UAV的位置和歐拉角信息，并且得到速度和加速度數(shù)據(jù)。在噪聲環(huán)境下，改進(jìn)差分進(jìn)化算法能很好地進(jìn)行參數(shù)中心點(diǎn)估計(jì)和優(yōu)化回歸函數(shù)參數(shù)。在集員區(qū)間估計(jì)方面，參數(shù) a1，a4，a5，a7，a9的估計(jì)區(qū)間較小，其估計(jì)結(jié)果可信度高，而剩下的5個(gè)參數(shù)估計(jì)區(qū)間較大，其結(jié)果只可作為參考，因?yàn)檫@5個(gè)參數(shù)數(shù)量級(jí)較小，即使在參數(shù)變化較大的情況下引起的誤差仍然較小，依然滿足可行解方程，所以區(qū)間估計(jì)較大。

6 結(jié)束語

本文僅基于UAV位置和歐拉角信息，利用跟蹤微分器得到速度和加速度數(shù)據(jù)，將復(fù)雜系統(tǒng)分解為兩個(gè)獨(dú)立的子系統(tǒng)進(jìn)行辨識(shí)，減輕了辨識(shí)的工作量。利用非線性向量回歸方法建立逼近參數(shù)向量與誤差向量的范數(shù)之間的函數(shù)關(guān)系模型，用差分進(jìn)化算法進(jìn)行參數(shù)中心點(diǎn)估計(jì)和優(yōu)化回歸函數(shù)參數(shù)，避免參數(shù)選擇的盲目性。為了進(jìn)一步完善算法，對(duì)樣本數(shù)據(jù)進(jìn)行加權(quán)處理。給出了一種四旋翼參數(shù)的集員估計(jì)方法，為集員辨識(shí)在飛行器上的應(yīng)用提供了一種新的思路。

［1］ Milanese M，Vicino A.Optimal estimation theory for dynamic system with set membership uncertainty:an overview ［J］.Automatica，1991，27(6):997-1009.

［2］ Cheng M F，Yurkovich S，Passino K M.An optimal volume ellipsoid algorithm for parameter set estimation［J］.IEEE Transactions on Automatic Control，1993，38(8):1292-1296.

［3］ Kosut R L，Lau M，Boyd S.Set-membership identification of systems with parametric and nonparametric uncertainty［J］.IEEE Transactions on Automatic Control，1992，37(7):929-941.

［4］ Fogel E，Huang Y F.On the value of information in system identificationˉbounded noise case［J］.Automatica，1982，18(2):229-238.

［5］ Deller J P.Set membership identification in digital signal processing［J］.IEEE ASSP Magazine，1989，6(10):4-20.

［6］ Veres SM，Norton JP.Structure identification of parameter-bounding models by use of noise－structure bounds［J］.International Journal of Control，1989，50(2):639-649.

［7］ Broman V，Shensa M J.A compact algorithm for the intersection and approximation of N-dimensional polytopes［J］.Mathematics and Computers in Simulation，1990，32(5):469-480.

［8］ Chakraborty U K，Das S，Konar A.Differential evolution with local neighborhood［C］//2006 IEEE Congress on Evolutionary Computation.Piscataway，N.J.IEEE，2006:2042-2049.

［9］ Wang X H，Chen Z Q，Yuan Z Z.Design and analysis for new discrete tracking-differentiators［J］.Applied MathematicˉA Journal of Chinese Universities，2003，18(2):214-222.

［10］ Suykens J A K，Vandewalle J.Least squares support vector machine classifiers ［J］.Neural Processing Letters，1999，9(3):293-300.

［11］ Chang W D.Parameter identification of Chen and Lü systems:a differential evolution approach ［J］.Chaos，Solitons ＆ Fractals，2007，32(4):1469-1476.

［12］ Das S，Konar A，Chakraborty U K.Two improved differential evolution schemes for faster global search［C］//Proceedings of the 2005 Conference on Genetic and Evolutionary Computation.New York，USA，ACM，2005:991-998.

［13］ Bingül Z，Karahan O.Dynamic identification of Staubli RX-60 robot using PSO and LS method［J］.Expert Systems with Applications，2011，38(4):4136-4149.

［14］ Tang H H，Xue ST，F(xiàn)an CX.Differential evolution strategy for structural system identification ［J］.Computers and Structures，2008，86(21):2004-2012.

［15］ Keesman K J，Stappers R J J.Nonlinear set-membership estimation:a support vector machine approach ［J］.Journal of Inverse and Ill-Posed Problems，2004，12(1):27-41.

［16］柴偉.集員估計(jì)理論、方法及其應(yīng)用［D］.北京:北京航空航天大學(xué)，2010.