一種基于模態參數實時辨識方法的參數時變航天器控制方法*

2019-12-03 12:09:50賈貴鵬趙育善

飛控與探測 2019年6期

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

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

0 引言

對于存在剛柔耦合航天器的姿態運動，為實現對姿態的控制，通常需要對系統狀態量進行估計，獲取系統盡可能精確的數學模型，但多數情況下，系統的數學精確模型不能通過理論計算獲得，甚至系統的參數存在時變，系統參數辨識就是確定系統數學模型的理論方法[1-2]。通過實驗分析進行參數辨識即通過選取結構上測試點，采集測試點的動態數據，通過特定的估計和辨識方法來獲取系統模型及參數。在辨識獲取系統姿態和模態狀態量后，即可利用無模型控制方法對航天器姿態進行控制，而不需要受控對象的其他模型參數。

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

本文首先建立參數時變剛柔耦合航天器的姿態動力學模型，其次利用文獻[11]的參數遞推辨識方法估算系統狀態參數，最后利用無模型控制方法對該系統進行控制，在不需要獲得系統其他模型參數的情況下可實現對系統的控制。

1 數學模型的建立

1.1 剛柔耦合航天器姿態時變動力學方程

對于帶有大型柔性帆板或柔性天線的航天器的建模問題，常用的方式為將其簡化為一個中心剛體附加若干柔性附件的模型結構，假設航天器帶有兩個可轉動的太陽能帆板和兩個具有桁架結構的天線，如圖1所示。

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

航天器本體為剛體，附件1、2為桁架結構的天線，附件3、4為太陽能帆板，帆板旋轉角標記為δ，當帆板面向xb軸方向時δ=0。假設帆板以很小的角速度旋轉以保持供能，可得到簡化的動力學方程[12]：

(1)

(2)

(3)

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

由于本文只研究剛柔耦合航天器的姿態控制，因此認為力F很小，航天器的線性運動對振動的影響也可以忽略不計，忽略航天器的線性運動，由式(2)、(3)可得以下航天器動力學方程：

(4)

(5)

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

簡化上述動力學方程，得到：

(6)

(7)

其中，

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

將上述動力學方程式寫成狀態空間方程形式：

(8)

(9)

(10)

量測方程可以寫為：

y(t)=Cx(t)

(11)

(12)

其中，C∈Rm×n為輸出矩陣；Φ為相應的模態矩陣。

2 遞推模態參數辨識方法

2.1 系統描述和假設

動態系統進行離散化并考慮噪聲對系統的影響，可以進一步將其寫為如下的新息形式：

xk+1=Akxk+Bkuk+Kkek

(13)

yk=Ckxk+ek

(14)

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

(15)

yk=Ckxk+ek

(16)

定義一個如下形式的VARX預測器，則k時刻的輸出yk可以表述為

(17)

yk=Ξkφk

(18)

式中：

那么利用自適應濾波技術，可以得到各個采樣時刻的矩陣Ξk的最小二乘遞推形式為：

(19)

(20)

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

2.2 基于仿射投影算法的系統狀態量估計方法

對于式(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算法，從而得到各個時刻的狀態量xk的遞推形式為：

(28)

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

(29)

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

(30)

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

(31)

(32)

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

(33)

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

3 基于辨識狀態量的無模型控制方法

無模型控制方法是一種數據驅動控制方法，根據辨識得到的狀態量xk可以利用無模型控制器對航天器進行控制，而不需要其他系統模型。考慮如下MIMO離散時間非線性系統

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離散時間非線性系統，提出如下兩個假設。

假設 3.1除有限時刻點外，fi(·),i=1,2，…,m，對各個變量的分量都存在連續偏導數。

(36)

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

考慮如下輸入準則函數

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

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

(37)

其中，λ>0是一個權重因子。

將式(36)代入準則函數(37)中，對u(k)求差分并令其等于零，得

Δu(k-i+1))

(38)

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

(39)

(40)

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

辨識與控制流程為：

(1)初始化

(2)辨識計算系統狀態量xk

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

(6)回到第(2)步進行下一步循環

4 仿真驗證

本文以文獻[11]中航天器的模態參數為初始值進行仿真實驗，本體初始姿態角為[0.2 0.3 0.3]rad，角速度為[0 0 0](rad/s)；太陽能帆板旋轉角速度為7.278×10-5(rad/s)(即360°/天)，四個附件分別選取的模態階數為8、8、7、7。對姿態角和四個附件的一階模態坐標的辨識結果如圖2所示。