一種基于模態(tài)參數(shù)實(shí)時(shí)辨識(shí)方法的參數(shù)時(shí)變航天器控制方法*

賈貴鵬， 趙 欣， 趙育善， 師 鵬

(北京航空航天大學(xué) 宇航學(xué)院·北京·100191)

0 引 言

對(duì)于存在剛?cè)狁詈虾教炱鞯淖藨B(tài)運(yùn)動(dòng)，為實(shí)現(xiàn)對(duì)姿態(tài)的控制，通常需要對(duì)系統(tǒng)狀態(tài)量進(jìn)行估計(jì)，獲取系統(tǒng)盡可能精確的數(shù)學(xué)模型，但多數(shù)情況下，系統(tǒng)的數(shù)學(xué)精確模型不能通過(guò)理論計(jì)算獲得，甚至系統(tǒng)的參數(shù)存在時(shí)變，系統(tǒng)參數(shù)辨識(shí)就是確定系統(tǒng)數(shù)學(xué)模型的理論方法[1-2]。通過(guò)實(shí)驗(yàn)分析進(jìn)行參數(shù)辨識(shí)即通過(guò)選取結(jié)構(gòu)上測(cè)試點(diǎn)，采集測(cè)試點(diǎn)的動(dòng)態(tài)數(shù)據(jù)，通過(guò)特定的估計(jì)和辨識(shí)方法來(lái)獲取系統(tǒng)模型及參數(shù)。在辨識(shí)獲取系統(tǒng)姿態(tài)和模態(tài)狀態(tài)量后，即可利用無(wú)模型控制方法對(duì)航天器姿態(tài)進(jìn)行控制，而不需要受控對(duì)象的其他模型參數(shù)。

在航天領(lǐng)域，為解決大型撓性航天器具有低頻密集的特點(diǎn)，時(shí)域辨識(shí)成為航天器參數(shù)辨識(shí)的主要方法，時(shí)域法無(wú)需激勵(lì)數(shù)據(jù)，直接利用響應(yīng)數(shù)據(jù)辨識(shí)系統(tǒng)參數(shù)，適合對(duì)系統(tǒng)進(jìn)行在線參數(shù)辨識(shí)，獲得系統(tǒng)動(dòng)態(tài)性能[3-7]。 Liu[8]利用基于奇異值分解(Singular Value Decomposition，SVD)方法獲取時(shí)變系統(tǒng)的參數(shù)估計(jì)，并提出偽模態(tài)方法來(lái)解決時(shí)變系統(tǒng)模態(tài)分析問(wèn)題。Shao等[9]將仿射投影算法(Affine Projection Algorithm，APA)和符號(hào)算法(Sign Algorithms，SAS)的優(yōu)點(diǎn)結(jié)合，通過(guò)使用多個(gè)投影同時(shí)根據(jù)L1范數(shù)優(yōu)化標(biāo)準(zhǔn)更新其權(quán)重，提出了仿射投影符號(hào)算法(Affine Projection Sign Algorithm，APSA)，該方法計(jì)算復(fù)雜度較低，抗脈沖干擾性強(qiáng)，收斂速度快。Houtzager等[10]提出基于遞歸預(yù)測(cè)器的子空間辨識(shí)法(Recursive Predictor-Based Subspace Identification，RPBSID)，通過(guò)自適應(yīng)濾波，利用外生向量自回歸(Vector Autoregressive with Exogenous，VARX)預(yù)測(cè)器對(duì)系統(tǒng)參數(shù)進(jìn)行辨識(shí)。針對(duì)RPBSID方法在對(duì)狀態(tài)量進(jìn)行估計(jì)時(shí)需要構(gòu)建廣義Hankel矩陣而導(dǎo)致計(jì)算量較大的問(wèn)題，倪智宇等[11]將RPBSID方法和APSA方法進(jìn)行結(jié)合，利用APSA方法對(duì)系統(tǒng)狀態(tài)量進(jìn)行估計(jì)，提高了計(jì)算效率。無(wú)模型控制是侯忠生提出的一種控制方法[12-14]，其主要特點(diǎn)是不需要系統(tǒng)的動(dòng)力學(xué)模型，只需要系統(tǒng)的激勵(lì)和響應(yīng)量即可對(duì)系統(tǒng)實(shí)現(xiàn)控制，上述研究對(duì)象集中在機(jī)械、化工等工業(yè)領(lǐng)域，對(duì)于復(fù)雜航天器控制問(wèn)題研究較少，本文針對(duì)一種因撓性結(jié)構(gòu)轉(zhuǎn)動(dòng)引起模態(tài)參數(shù)變化的航天器的模態(tài)辨識(shí)與控制問(wèn)題進(jìn)行了研究。

本文首先建立參數(shù)時(shí)變剛?cè)狁詈虾教炱鞯淖藨B(tài)動(dòng)力學(xué)模型，其次利用文獻(xiàn)[11]的參數(shù)遞推辨識(shí)方法估算系統(tǒng)狀態(tài)參數(shù)，最后利用無(wú)模型控制方法對(duì)該系統(tǒng)進(jìn)行控制，在不需要獲得系統(tǒng)其他模型參數(shù)的情況下可實(shí)現(xiàn)對(duì)系統(tǒng)的控制。

1 數(shù)學(xué)模型的建立

1.1 剛?cè)狁詈虾教炱髯藨B(tài)時(shí)變動(dòng)力學(xué)方程

對(duì)于帶有大型柔性帆板或柔性天線的航天器的建模問(wèn)題，常用的方式為將其簡(jiǎn)化為一個(gè)中心剛體附加若干柔性附件的模型結(jié)構(gòu)，假設(shè)航天器帶有兩個(gè)可轉(zhuǎn)動(dòng)的太陽(yáng)能帆板和兩個(gè)具有桁架結(jié)構(gòu)的天線，如圖1所示。

圖1 航天器本體和附件坐標(biāo)系示意圖Fig.1 Body coordinate system of the spacecraft

航天器本體為剛體，附件1、2為桁架結(jié)構(gòu)的天線，附件3、4為太陽(yáng)能帆板，帆板旋轉(zhuǎn)角標(biāo)記為δ，當(dāng)帆板面向xb軸方向時(shí)δ=0。假設(shè)帆板以很小的角速度旋轉(zhuǎn)以保持供能，可得到簡(jiǎn)化的動(dòng)力學(xué)方程[12]：

(1)

(2)

(3)

式中，r=[rx,ry,rz]T為 J2000慣性坐標(biāo)系下系統(tǒng)質(zhì)心的位置向量，θ=[θroll,θpitch,θyaw]為姿態(tài)角矢量，ηi∈Rni為第i個(gè)撓性附件的模態(tài)坐標(biāo)，Zi=diag[ζai,1,…,ζai,2]∈Rni*ni為對(duì)應(yīng)ηi的阻尼比矩陣，Ωi=diag[ωai,1,…,ωai,2]∈Rni*ni為對(duì)應(yīng)ηi的模態(tài)頻率矩陣，F(xiàn)∈R3為作用在航天器系統(tǒng)上的力，T∈R3為作用在航天器系統(tǒng)上的力矩，M,J,Υ,Γi,Pi的定義如下：

由于本文只研究剛?cè)狁詈虾教炱鞯淖藨B(tài)控制，因此認(rèn)為力F很小，航天器的線性運(yùn)動(dòng)對(duì)振動(dòng)的影響也可以忽略不計(jì)，忽略航天器的線性運(yùn)動(dòng)，由式(2)、(3)可得以下航天器動(dòng)力學(xué)方程：

(4)

(5)

其中，Λi=2ZiΩi。假設(shè)帆板角度δ是時(shí)變的，即δ(t)，且J和P的變化依賴于帆板轉(zhuǎn)角δ(t)。故可用J(t)和Pi(t)取代J和P。T為由反作用飛輪和擾動(dòng)輸入產(chǎn)生的三軸控制力矩輸入。

簡(jiǎn)化上述動(dòng)力學(xué)方程，得到：

(6)

(7)

其中，

其中，u∈R3為三軸控制力矩輸入，w∈R3為干擾力矩輸入。

將上述動(dòng)力學(xué)方程式寫(xiě)成狀態(tài)空間方程形式：

(8)

(9)

(10)

量測(cè)方程可以寫(xiě)為：

y(t)=Cx(t)

(11)

(12)

其中，C∈Rm*n為輸出矩陣；Φ為相應(yīng)的模態(tài)矩陣。

2 遞推模態(tài)參數(shù)辨識(shí)方法

2.1 系統(tǒng)描述和假設(shè)

動(dòng)態(tài)系統(tǒng)進(jìn)行離散化并考慮噪聲對(duì)系統(tǒng)的影響，可以進(jìn)一步將其寫(xiě)為如下的新息形式：

xk+1=Akxk+Bkuk+Kkek

(13)

yk=Ckxk+ek

(14)

其中，xk∈Rn*1為狀態(tài)向量；uk∈Rr*1為輸入向量；yk∈Rm*1為輸出向量；下標(biāo)k為離散后的時(shí)間點(diǎn)，ek表示新息白噪聲序列，而矩陣Kk表示卡爾曼增益矩陣；A∈Rn*n為狀態(tài)矩陣；B∈Rn*r為輸入矩陣；C∈Rm*n為輸出矩陣。 那么對(duì)式(13)和(14)做進(jìn)一步的變換，可以寫(xiě)為：

(15)

yk=Ckxk+ek

(16)

定義一個(gè)如下形式的VARX預(yù)測(cè)器，則k時(shí)刻的輸出yk可以表述為

(17)

yk=Ξkφk

(18)

式中：

那么利用自適應(yīng)濾波技術(shù)，可以得到各個(gè)采樣時(shí)刻的矩陣Ξk的最小二乘遞推形式為：

(19)

(20)

式中：β1為遺忘因子，需要保證0<β1≤1；Ζk為定義的臨時(shí)變量矩陣，通常初值給為Z0=I。通過(guò)式(19)的遞推計(jì)算，可以得到各個(gè)時(shí)刻的Markov參數(shù)矩陣Ξk的值，接下來(lái)利用該矩陣遞推估計(jì)系統(tǒng)的狀態(tài)量xk。

2.2 基于仿射投影算法的系統(tǒng)狀態(tài)量估計(jì)方法

對(duì)于式(15)、(16)，如果暫不考慮量測(cè)噪聲ek，那么從k-p到k時(shí)刻的輸出信號(hào){yk-p,yk-p+1,…,yk}可以分別寫(xiě)為：

yk-p=Ck-pxk-p

(21)

(22)

(23)

則式(21)～(23)可以用矩陣形式表示為：

(24)

式中：

τk-p=Γk-pxk-p

(25)

那么對(duì)于k-1時(shí)刻，式(25)可以表示為如下的最小二乘形式：

τk-1=Γk-1xk-1

(26)

其中

(27)

對(duì)于式(26)，這里引入APA算法，從而得到各個(gè)時(shí)刻的狀態(tài)量xk的遞推形式為：

(28)

Φk-1=τk-1-Γk-1xk-1

(29)

式中：μ為遞推過(guò)程中的仿射投影因子，一般應(yīng)滿足0<μ≤1。在式(28)的遞推迭代中，矩陣Γk-1的初值Γ0應(yīng)保證為滿秩矩陣，如果沒(méi)有其他先驗(yàn)知識(shí)，則Γ0一般可以選擇為

(30)

在式(28)中，在遞推得到xk后，矩陣Γk的更新同樣可以通過(guò)如下的遞推最小二乘算法得到：

(31)

(32)

Γk=Γk-1+(τk-Γk-1xk)Wk

(33)

式中：β2為遺忘因子，而矩陣Lk的初值選取方式和Γ0相似。根據(jù)式(28)～(33)，可以遞推計(jì)算得到各個(gè)時(shí)刻的狀態(tài)量xk。

3 基于辨識(shí)狀態(tài)量的無(wú)模型控制方法

無(wú)模型控制方法是一種數(shù)據(jù)驅(qū)動(dòng)控制方法，根據(jù)辨識(shí)得到的狀態(tài)量xk可以利用無(wú)模型控制器對(duì)航天器進(jìn)行控制，而不需要其他系統(tǒng)模型?？紤]如下MIMO離散時(shí)間非線性系統(tǒng)

y(k+1)=

f(y(k),…,y(k-ny),u(k),…,u(k-nu))

(34)

uT(k),…,uT(k-Lu+1)]T

(35)

對(duì)形如式(34)的MIMO離散時(shí)間非線性系統(tǒng)，提出如下兩個(gè)假設(shè)。

假設(shè) 3.1除有限時(shí)刻點(diǎn)外，fi(·),i=1,2，…,m，對(duì)各個(gè)變量的分量都存在連續(xù)偏導(dǎo)數(shù)。

(36)

且對(duì)任意時(shí)刻k，Φf,Ly,Lu(k)=[Φ1(k)…ΦLy+Lu(k)] 是有界的，其中Φi(k)∈Rm*m,i=1,…Ly+Lu。

考慮如下輸入準(zhǔn)則函數(shù)

J(u(k))=‖y*(k+1)-y(k+1)‖2+

λ‖u(k)-u(k-1)‖2

(37)

其中，λ>0是一個(gè)權(quán)重因子。

將式(36)代入準(zhǔn)則函數(shù)(37)中，對(duì)u(k)求差分并令其等于零，得

Δu(k-i+1))

(38)

考慮到該控制算法中包含了矩陣求逆運(yùn)算，在矩陣維數(shù)較大時(shí)運(yùn)算量大，運(yùn)算速度較慢，根據(jù)無(wú)模型自適應(yīng)控制(Model Free Adaptive Control，MFAC)理論，可以得到針對(duì)離散時(shí)間MIMO非線性系統(tǒng)的MFAC控制方案如下

(39)

(40)

其中，i=1,2，…,m;j=1,2，…,m;ρ1,ρ2,…ρLy+Lu∈(0,1]；η∈(0,2]；λ>0，μ>0。

辨識(shí)與控制流程為：

(1)初始化

(2)辨識(shí)計(jì)算系統(tǒng)狀態(tài)量xk

(5)由(40)計(jì)算控制輸入u(k)

(6)回到第(2)步進(jìn)行下一步循環(huán)

4 仿真驗(yàn)證

本文以文獻(xiàn)[11]中航天器的模態(tài)參數(shù)為初始值進(jìn)行仿真實(shí)驗(yàn)，本體初始姿態(tài)角為[0.2 0.3 0.3]rad，角速度為[0 0 0](rad/s)；太陽(yáng)能帆板旋轉(zhuǎn)角速度為7.278*10-5(rad/s)(即360°/天)，四個(gè)附件分別選取的模態(tài)階數(shù)為8、8、7、7。對(duì)姿態(tài)角和四個(gè)附件的一階模態(tài)坐標(biāo)的辨識(shí)結(jié)果如圖2所示。

(a)姿態(tài)角φ辨識(shí)結(jié)果

(b)姿態(tài)角θ辨識(shí)結(jié)果

(c)姿態(tài)角ψ辨識(shí)結(jié)果

(d)附件1的1階模態(tài)坐標(biāo)辨識(shí)結(jié)果

(e)附件2的1階模態(tài)坐標(biāo)辨識(shí)結(jié)果

(f)附件3的1階模態(tài)坐標(biāo)辨識(shí)結(jié)果

(g)附件4的1階模態(tài)坐標(biāo)辨識(shí)結(jié)果圖2 歐拉角與四個(gè)附件的1階模態(tài)坐標(biāo)的辨識(shí)結(jié)果Fig.2 Identification result of Euler angles and firstorder modal coordinate of four appendage

根據(jù)實(shí)驗(yàn)結(jié)果可以看出，仿真初始階段的辨識(shí)誤差較大，隨著辨識(shí)過(guò)程的繼續(xù)，辨識(shí)誤差趨于減小。在參數(shù)發(fā)生劇烈變化時(shí)，辨識(shí)誤差相對(duì)較大。

利用設(shè)計(jì)的辨識(shí)與控制方案對(duì)上述航天器進(jìn)行控制，設(shè)姿態(tài)角為[0.2 0.3 0.3]rad，角速度為[0.01 -0.01 0.01](rad/s)，控制器參數(shù)取ρ=[0.5 0.5]，Lu=Ly=1，η=0.1，λ=1，μ=1。仿真結(jié)果如圖3所示。

(a)姿態(tài)角控制結(jié)果

(b)角速度控制結(jié)果

(c)控制力矩仿真結(jié)果圖3 MFAC控制效果圖Fig. 3 MFAC control effect

觀察仿真結(jié)果，初始階段由于辨識(shí)誤差較大，控制過(guò)程的波動(dòng)較大，隨著辨識(shí)過(guò)程的繼續(xù)，航天器姿態(tài)最終趨于期望值。

4 結(jié) 論

本文研究了剛?cè)狁詈虾教炱鳡顟B(tài)估計(jì)的問(wèn)題。本文采用改進(jìn)RPBSID來(lái)實(shí)現(xiàn)狀態(tài)量的遞推辨識(shí)，該方法只需要利用系統(tǒng)的激勵(lì)響應(yīng)信息來(lái)估計(jì)系統(tǒng)狀態(tài)量，利用估計(jì)的狀態(tài)量采用無(wú)模型控制方法對(duì)系統(tǒng)進(jìn)行控制，在不需要系統(tǒng)真實(shí)動(dòng)力學(xué)模型參數(shù)的情況下可實(shí)現(xiàn)對(duì)系統(tǒng)的控制。通過(guò)仿真實(shí)驗(yàn)，驗(yàn)證了辨識(shí)方法和控制方法的有效性。