賈貴鵬, 趙 欣, 趙育善, 師 鵬
(北京航空航天大學(xué) 宇航學(xué)院·北京·100191)
對于存在剛?cè)狁詈虾教炱鞯淖藨B(tài)運動,為實現(xiàn)對姿態(tài)的控制,通常需要對系統(tǒng)狀態(tài)量進行估計,獲取系統(tǒng)盡可能精確的數(shù)學(xué)模型,但多數(shù)情況下,系統(tǒng)的數(shù)學(xué)精確模型不能通過理論計算獲得,甚至系統(tǒng)的參數(shù)存在時變,系統(tǒng)參數(shù)辨識就是確定系統(tǒng)數(shù)學(xué)模型的理論方法[1-2]。通過實驗分析進行參數(shù)辨識即通過選取結(jié)構(gòu)上測試點,采集測試點的動態(tài)數(shù)據(jù),通過特定的估計和辨識方法來獲取系統(tǒng)模型及參數(shù)。在辨識獲取系統(tǒng)姿態(tài)和模態(tài)狀態(tài)量后,即可利用無模型控制方法對航天器姿態(tài)進行控制,而不需要受控對象的其他模型參數(shù)。
在航天領(lǐng)域,為解決大型撓性航天器具有低頻密集的特點,時域辨識成為航天器參數(shù)辨識的主要方法,時域法無需激勵數(shù)據(jù),直接利用響應(yīng)數(shù)據(jù)辨識系統(tǒng)參數(shù),適合對系統(tǒng)進行在線參數(shù)辨識,獲得系統(tǒng)動態(tài)性能[3-7]。 Liu[8]利用基于奇異值分解(Singular Value Decomposition,SVD)方法獲取時變系統(tǒng)的參數(shù)估計,并提出偽模態(tài)方法來解決時變系統(tǒng)模態(tài)分析問題。Shao等[9]將仿射投影算法(Affine Projection Algorithm,APA)和符號算法(Sign Algorithms,SAS)的優(yōu)點結(jié)合,通過使用多個投影同時根據(jù)L1范數(shù)優(yōu)化標準更新其權(quán)重,提出了仿射投影符號算法(Affine Projection Sign Algorithm,APSA),該方法計算復(fù)雜度較低,抗脈沖干擾性強,收斂速度快。Houtzager等[10]提出基于遞歸預(yù)測器的子空間辨識法(Recursive Predictor-Based Subspace Identification,RPBSID),通過自適應(yīng)濾波,利用外生向量自回歸(Vector Autoregressive with Exogenous,VARX)預(yù)測器對系統(tǒng)參數(shù)進行辨識。針對RPBSID方法在對狀態(tài)量進行估計時需要構(gòu)建廣義Hankel矩陣而導(dǎo)致計算量較大的問題,倪智宇等[11]將RPBSID方法和APSA方法進行結(jié)合,利用APSA方法對系統(tǒng)狀態(tài)量進行估計,提高了計算效率。無模型控制是侯忠生提出的一種控制方法[12-14],其主要特點是不需要系統(tǒng)的動力學(xué)模型,只需要系統(tǒng)的激勵和響應(yīng)量即可對系統(tǒng)實現(xiàn)控制,上述研究對象集中在機械、化工等工業(yè)領(lǐng)域,對于復(fù)雜航天器控制問題研究較少,本文針對一種因撓性結(jié)構(gòu)轉(zhuǎn)動引起模態(tài)參數(shù)變化的航天器的模態(tài)辨識與控制問題進行了研究。
本文首先建立參數(shù)時變剛?cè)狁詈虾教炱鞯淖藨B(tài)動力學(xué)模型,其次利用文獻[11]的參數(shù)遞推辨識方法估算系統(tǒng)狀態(tài)參數(shù),最后利用無模型控制方法對該系統(tǒng)進行控制,在不需要獲得系統(tǒng)其他模型參數(shù)的情況下可實現(xiàn)對系統(tǒng)的控制。
對于帶有大型柔性帆板或柔性天線的航天器的建模問題,常用的方式為將其簡化為一個中心剛體附加若干柔性附件的模型結(jié)構(gòu),假設(shè)航天器帶有兩個可轉(zhuǎn)動的太陽能帆板和兩個具有桁架結(jié)構(gòu)的天線,如圖1所示。
圖1 航天器本體和附件坐標系示意圖Fig.1 Body coordinate system of the spacecraft
航天器本體為剛體,附件1、2為桁架結(jié)構(gòu)的天線,附件3、4為太陽能帆板,帆板旋轉(zhuǎn)角標記為δ,當帆板面向xb軸方向時δ=0。假設(shè)帆板以很小的角速度旋轉(zhuǎn)以保持供能,可得到簡化的動力學(xué)方程[12]:
(1)
(2)
(3)
式中,r=[rx,ry,rz]T為 J2000慣性坐標系下系統(tǒng)質(zhì)心的位置向量,θ=[θroll,θpitch,θyaw]為姿態(tài)角矢量,ηi∈Rni為第i個撓性附件的模態(tài)坐標,Zi=diag[ζai,1,…,ζai,2]∈Rni×ni為對應(yīng)ηi的阻尼比矩陣,Ωi=diag[ωai,1,…,ωai,2]∈Rni×ni為對應(yīng)ηi的模態(tài)頻率矩陣,F(xiàn)∈R3為作用在航天器系統(tǒng)上的力,T∈R3為作用在航天器系統(tǒng)上的力矩,M,J,Υ,Γi,Pi的定義如下:
由于本文只研究剛?cè)狁詈虾教炱鞯淖藨B(tài)控制,因此認為力F很小,航天器的線性運動對振動的影響也可以忽略不計,忽略航天器的線性運動,由式(2)、(3)可得以下航天器動力學(xué)方程:
(4)
(5)
其中,Λi=2ZiΩi。假設(shè)帆板角度δ是時變的,即δ(t),且J和P的變化依賴于帆板轉(zhuǎn)角δ(t)。故可用J(t)和Pi(t)取代J和P。T為由反作用飛輪和擾動輸入產(chǎn)生的三軸控制力矩輸入。
簡化上述動力學(xué)方程,得到:
(6)
(7)
其中,
其中,u∈R3為三軸控制力矩輸入,w∈R3為干擾力矩輸入。
將上述動力學(xué)方程式寫成狀態(tài)空間方程形式:
(8)
(9)
(10)
量測方程可以寫為:
y(t)=Cx(t)
(11)
(12)
其中,C∈Rm×n為輸出矩陣;Φ為相應(yīng)的模態(tài)矩陣。
動態(tài)系統(tǒng)進行離散化并考慮噪聲對系統(tǒng)的影響,可以進一步將其寫為如下的新息形式:
xk+1=Akxk+Bkuk+Kkek
(13)
yk=Ckxk+ek
(14)
其中,xk∈Rn×1為狀態(tài)向量;uk∈Rr×1為輸入向量;yk∈Rm×1為輸出向量;下標k為離散后的時間點,ek表示新息白噪聲序列,而矩陣Kk表示卡爾曼增益矩陣;A∈Rn×n為狀態(tài)矩陣;B∈Rn×r為輸入矩陣;C∈Rm×n為輸出矩陣。 那么對式(13)和(14)做進一步的變換,可以寫為:
(15)
yk=Ckxk+ek
(16)
定義一個如下形式的VARX預(yù)測器,則k時刻的輸出yk可以表述為
(17)
yk=Ξkφk
(18)
式中:
那么利用自適應(yīng)濾波技術(shù),可以得到各個采樣時刻的矩陣Ξk的最小二乘遞推形式為:
(19)
(20)
式中:β1為遺忘因子,需要保證0<β1≤1;Ζk為定義的臨時變量矩陣,通常初值給為Z0=I。通過式(19)的遞推計算,可以得到各個時刻的Markov參數(shù)矩陣Ξk的值,接下來利用該矩陣遞推估計系統(tǒng)的狀態(tài)量xk。
對于式(15)、(16),如果暫不考慮量測噪聲ek,那么從k-p到k時刻的輸出信號{yk-p,yk-p+1,…,yk}可以分別寫為:
yk-p=Ck-pxk-p
(21)
(22)
(23)
則式(21)~(23)可以用矩陣形式表示為:
(24)
式中:
τk-p=Γk-pxk-p
(25)
那么對于k-1時刻,式(25)可以表示為如下的最小二乘形式:
τk-1=Γk-1xk-1
(26)
其中
(27)
對于式(26),這里引入APA算法,從而得到各個時刻的狀態(tài)量xk的遞推形式為:
(28)
Φk-1=τk-1-Γk-1xk-1
(29)
式中:μ為遞推過程中的仿射投影因子,一般應(yīng)滿足0<μ≤1。在式(28)的遞推迭代中,矩陣Γk-1的初值Γ0應(yīng)保證為滿秩矩陣,如果沒有其他先驗知識,則Γ0一般可以選擇為
(30)
在式(28)中,在遞推得到xk后,矩陣Γk的更新同樣可以通過如下的遞推最小二乘算法得到:
(31)
(32)
Γk=Γk-1+(τk-Γk-1xk)Wk
(33)
式中:β2為遺忘因子,而矩陣Lk的初值選取方式和Γ0相似。根據(jù)式(28)~(33),可以遞推計算得到各個時刻的狀態(tài)量xk。
無模型控制方法是一種數(shù)據(jù)驅(qū)動控制方法,根據(jù)辨識得到的狀態(tài)量xk可以利用無模型控制器對航天器進行控制,而不需要其他系統(tǒng)模型。考慮如下MIMO離散時間非線性系統(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)
對形如式(34)的MIMO離散時間非線性系統(tǒng),提出如下兩個假設(shè)。
假設(shè) 3.1除有限時刻點外,fi(·),i=1,2,…,m,對各個變量的分量都存在連續(xù)偏導(dǎo)數(shù)。
(36)
且對任意時刻k,Φf,Ly,Lu(k)=[Φ1(k)…ΦLy+Lu(k)] 是有界的,其中Φi(k)∈Rm×m,i=1,…Ly+Lu。
考慮如下輸入準則函數(shù)
J(u(k))=||y*(k+1)-y(k+1)||2+
λ||u(k)-u(k-1)||2
(37)
其中,λ>0是一個權(quán)重因子。
將式(36)代入準則函數(shù)(37)中,對u(k)求差分并令其等于零,得
Δu(k-i+1))
(38)
考慮到該控制算法中包含了矩陣求逆運算,在矩陣維數(shù)較大時運算量大,運算速度較慢,根據(jù)無模型自適應(yīng)控制(Model Free Adaptive Control,MFAC)理論,可以得到針對離散時間MIMO非線性系統(tǒng)的MFAC控制方案如下
(39)
(40)
其中,i=1,2,…,m;j=1,2,…,m;ρ1,ρ2,…ρLy+Lu∈(0,1];η∈(0,2];λ>0,μ>0。
辨識與控制流程為:
(1)初始化
(2)辨識計算系統(tǒng)狀態(tài)量xk
(5)由(40)計算控制輸入u(k)
(6)回到第(2)步進行下一步循環(huán)
本文以文獻[11]中航天器的模態(tài)參數(shù)為初始值進行仿真實驗,本體初始姿態(tài)角為[0.2 0.3 0.3]rad,角速度為[0 0 0](rad/s);太陽能帆板旋轉(zhuǎn)角速度為7.278×10-5(rad/s)(即360°/天),四個附件分別選取的模態(tài)階數(shù)為8、8、7、7。對姿態(tài)角和四個附件的一階模態(tài)坐標的辨識結(jié)果如圖2所示。
(a)姿態(tài)角φ辨識結(jié)果
(b)姿態(tài)角θ辨識結(jié)果
(c)姿態(tài)角ψ辨識結(jié)果
(d)附件1的1階模態(tài)坐標辨識結(jié)果
(e)附件2的1階模態(tài)坐標辨識結(jié)果
(f)附件3的1階模態(tài)坐標辨識結(jié)果
(g)附件4的1階模態(tài)坐標辨識結(jié)果圖2 歐拉角與四個附件的1階模態(tài)坐標的辨識結(jié)果Fig.2 Identification result of Euler angles and first- order modal coordinate of four appendage
根據(jù)實驗結(jié)果可以看出,仿真初始階段的辨識誤差較大,隨著辨識過程的繼續(xù),辨識誤差趨于減小。在參數(shù)發(fā)生劇烈變化時,辨識誤差相對較大。
利用設(shè)計的辨識與控制方案對上述航天器進行控制,設(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é)果,初始階段由于辨識誤差較大,控制過程的波動較大,隨著辨識過程的繼續(xù),航天器姿態(tài)最終趨于期望值。
本文研究了剛?cè)狁詈虾教炱鳡顟B(tài)估計的問題。本文采用改進RPBSID來實現(xiàn)狀態(tài)量的遞推辨識,該方法只需要利用系統(tǒng)的激勵響應(yīng)信息來估計系統(tǒng)狀態(tài)量,利用估計的狀態(tài)量采用無模型控制方法對系統(tǒng)進行控制,在不需要系統(tǒng)真實動力學(xué)模型參數(shù)的情況下可實現(xiàn)對系統(tǒng)的控制。通過仿真實驗,驗證了辨識方法和控制方法的有效性。