輪系航天器的角動量包絡分析及角動量管理

，

1. 北京控制工程研究所，北京 100094 2. 空間智能控制技術重點實驗室，北京 100094

反作用輪由于高精度連續方式輸出力矩，因此廣泛地作為航天器姿態控制系統的執行機構。為實現三軸姿態控制以及滿足角動量需求且具有一定冗余，一般配置4個或以上飛輪組成輪系系統輸出控制力矩[1-2]。航天器在軌運行時受環境外擾力矩作用，當無角動量卸載介入時輪系積累角動量將持續增長，使得飛輪組角動量飽和而喪失三軸控制能力。因此，系統一般配置磁卸載或噴氣卸載裝置[3]，特別是當判定飛輪角動量超過設定安全閾值時自主觸發噴氣卸載。然而，噴氣卸載勢必引起星體姿態波動對載荷工作產生影響，因此合理的系統設計，并輔之有效的角動量管理以減小不必要噴氣卸載，是在軌航天器實現載荷長時間連續工作的重要手段。

對于輪系航天器控制系統，輪系角動量包絡分析是輪系構型、飛輪容量選擇及角動量管理[3-4]，乃至輪系可重構性分析[5]等的先決條件。文獻[6]對飛輪角動量包絡研究首次闡述了實際最大角動量和偽逆解可達角動量兩類角動量包絡概念，利用幾何方法分析對比了幾種典型飛輪構型的兩類角動量包絡情況，并給出了一種斜裝飛輪系統基于偽逆解可達角動量包絡的卸載方案。文獻[7]針對同類型的三飛輪及四飛輪系統可能存在的不同構型，從最小能耗角度分析合適的安裝方式。

針對現常用能量最優的偽逆力矩分配方法未能充分利用輪系角動量和力矩的問題，相關研究[8-10]基于角動量與力矩包絡分析及新型優化力矩分配方法進行研究。Markley等基于幾何方法對各種構型的最大力矩及最大角動量包絡分析，并將包絡面分成與飛輪狀態相關的多區域，通過判定期望力矩所屬區域而采用預先按區域設定的不同L∞力矩分配陣進行力矩分配，以發揮輪系最大能力[8]。針對Markley給出包絡確定及解算效率，Yoon等采用包絡面上飛輪飽和狀態的角動量包絡及力矩包絡描述方式，結合包絡分析開展適應姿態機動的構型優化，并改進L∞力矩分配算法[9]。郭延寧等主要研究了航天器姿控系統中冗余飛輪構型的力矩分配問題，給出了能量最優的力矩分配策略,基于靜態最優化理論和輪系性質設計了相應的力矩最優分配策略,以有效地輸出飛輪系的最大力矩[10]。

角動量噴氣卸載除了需消耗有限的工質外，還可能對衛星軌道產生影響[11]，因此如何避免在特定力矩分配下系統角動量在達到包絡面之前出現飛輪角動量飽和而觸發不必要動量卸載的問題，一直是高軌航天器控制系統設計中需要格外關注的問題。針對既定任意冗余輪系構型航天器，本文基于可視化角動量包絡分析及合適角動量包絡特征參數描述選取研究，提出結合偽逆與基于目標角動量的零運動動態調整相結合的力矩分配方法，使得系統可充分利用實際最大角動量包絡，且適合航天器在軌實時有效的角動量管理。

1 輪系角動量包絡分析

考慮任意兩飛輪安裝不平行的N(N≥3)個飛輪所組成輪系，飛輪i(i=1,2,…,N)安裝方向hw,i∈R3為單位向量，對應角動量Hw,i∈R滿足如下約束：

式中：Hmin ,i與Hmax ,i為飛輪角動量最小值與最大值，當采用同類飛輪時一般有Hmax ,i=-Hmin ,i>0。

記角動量向量為：

輪系合成角動量可表示為:

對于星體坐標系下的任意單位向量

可使得S為：

記HΣ為輪系合成角動量HΣ在S方向的投影，則有

上式兩邊左乘矩陣C，有

式(6)可分解為如下兩部分：

求輪系組合沿星體S方向的最大合成角動量問題,即為在不等式(1)和等式(7)約束下的線性規劃問題:

以4個角動量相同的飛輪所組成的如圖1所示金字塔構型為例，飛輪角動量單位化且滿足-1≤Hw,i≤1(i=1,2,3，4)，各飛輪安裝方向為：

hw，1=[ 0.8165408118 0 0.5772877121 ]T

hw，2=[ 0 0.8165408118 0.5772877121 ]T

hw，3=[ -0.8165408118 0 0.5772877121 ]T

hw，4=[ 0 -0.8165408118 0.5772877121]T

對上述線性優化問題求解，可求得各飛輪角動量滿足Hmax ,i=-Hmin ,i=1 N·m·s條件下的輪系角動量包絡，如圖2所示。由可視化包絡圖可知輪系角動量包絡為十二面體，即包絡面由12個平面四邊形所圍成的封閉面。

綜合輪系合成角動量在可視化包絡面所對應的各飛輪角動量值可直觀地表現出角動量包絡的相關性質[8]：輪系合成角動量可看作由N維超立方體映射到三維空間的多面體，該多面體由頂點、邊和面組成，其中頂點為所有飛輪飽和角動量合成，各邊為除一個飛輪外其余飛輪均飽和下的角動量合成，面為除兩飛輪外其余飛輪均飽和下的角動量合成。

2 基于包絡面可視化的輪系角動 量分析

根據每個飛輪角動量約束條件式(1)的緊集特性可得Hw,i為凸集，進一步根據凸集相關性質[10]可得飛輪角動量包絡面的相關性質：輪系合成角動量HΣ為各個飛輪角動量Hw,i的線性組合，根據凸集之和為凸集的性質可知系統合成角動量HΣ形成的動量體也為凸集。另外，由于空間平面為凸集，于是根據凸集交集為凸集的性質可知角動量體外包絡平面多邊形為凸多邊形。

對于任意兩飛輪，不妨記為k與l，當其他N-2個飛輪角動量均處于飽和，Hw,k與Hw,l在其范圍內變化時所對應輪系合成角動量變化范圍為一個平行四邊形區域，其法線為n=hw,k×hw,l/‖hw,k×hw,l‖，且兩邊分別與hw,k和hw,l平行(文中‖·‖均為2-范數，表示向量長度)。當N-2個飽和飛輪的合成角動量在法線n上投影達到最大(或最小)時所對應平行四邊形則在輪系角動量包絡上，不妨記為Π:A1A2A3A4，如圖3所示，4個頂點按平面法線n(取為外法線)右手螺旋依次為A1～A4，Ai坐標為(xiyizi)。對于N個飛輪組成輪系，由任意兩飛輪組合可得包絡上平行四邊形數為N(N-1)。

由多邊形Π的任意第i個頂點Ai為起點，與其相鄰兩頂點Ai+1、Ai-1為終點，可形成兩個向量Ai，i+1與Ai，i-1(當i=4時為A4,1與A4,3)。由其中任意兩相鄰向量可求得平面Π的外法線n為：

結合任意一頂點Ai，可得Π所在平面的方程為：

x=Sxd，y=Syd，z=Szd(11)

式中：d為參量。

將直線方程各坐標系分量代入式(10)可得關于參量d的方程為：

nx(Sxd-xi)+ny(Syd-yi)+

且當nTS≠0時，其解為：

式中：AO，1為以原點O為起點、頂點A1為終點的向量。當d>0時意味著直線由原點沿單位向量S與四邊形Π所在平面存在交點，且向量B=d·S。

當判定沿單位向量S與平面四邊形Π所在平面存在相交點后，可進一步依據平面凸多邊形內點判斷方法[14]直接導出的沿方向S與包絡面上平面四邊形Π的交點存在性判據進行判斷。

交點存在性判據：若式(12)存在參量d大于零的解，且對于平面四邊形Π各邊均滿足如下不等式：

式中：Ai,B為以頂點Ai為起點、以交點B為終點的向量，i=1,2,3,4。則沿S方向與輪系角動量包絡面的四邊形Π存在交點B，當有等號成立時則交點B為Π邊上的點；否則與平面四邊形Π不相交。

由前可知，根據N飛輪組的輪系角動量包絡面由包絡多面體的其中4個頂點圍成的p=N(N-1)個平行四邊形Πυ(υ=1,2，…,p)區域所組成，并且當已知各平行四邊形Πυ內飛輪角動量飽和狀態即可求包絡面各點處對應飛輪極限角動量，因此采用如下3個包絡特征參量集的定義可對輪系角動量包絡進行描述：

1){Πυ:υ∈{1,2，…,p}}；

2)Πυ的有序頂點集{Λυ,j:j∈{1,2，3,4}}；

3) 包絡面上飛輪飽和狀態{Υυ,k:k∈{1,

2，…,N}}

注1：由于不同平面四邊形可能存在公共頂點，在具體實現上為減小存貯量可按一定編號存儲角動量包絡多面體頂點坐標，Πυ的頂點集各坐標可按索引方式獲??；當Hmax ,i=-Hmin ,i時，由式(2)有HΣ(-Hw)=-HΣ(Hw)，即包絡面上任一平行四邊形均存在與之原點對稱的平行四邊形，因而該情況下僅存貯一半數目平行四邊形的相關信息即可。

注2：為減少計算量，由可視化包絡可得合成角動量最大幅值HΣmax，當式(13)的解滿足0≤d≤HΣmax時才執行后續交點存在性判斷，否則直接判定該交點不在當前平行四邊形內。

以第1節的金字塔飛輪組合構型，由圖2所示可視化角動量包絡可得其12個面的包絡特征參數如表1所示。對于指定單位方向

按本節給出的過程可求得沿S方向與包絡面相交點在平行四邊形Π2內，且此處對應各飛輪角動量分別為1.00，1.00，-0.09和0.51，與由式(8)優化求解結果一致。

表1 金字塔構型的動量包絡面特征參數Table 1 Characteristic parameters of momentum envelope for RWs in pyramid array

3 輪系航天器的角動量管理

已知星體相對參考系的姿態誤差四元數Δq(t)和角速度誤差Δω(t)，類似文獻[15]可采用如下姿態控制律:

τc=KpΔq+KdΔω+ω×(Jsω+HΣ)

(15)

式中：Js∈R3×3為星體的轉動慣量；ω∈R3為星體角速度，ω×∈R3×3為ω的反對稱陣；控制系數陣Kp、Kd∈R3×3為正定矩陣。

假設N個飛輪組的當前角動量向量為Hw，且設定的飛輪偏置標稱角動量向量為Hw0∈RN。不失一般性，假設CwHw0=0，則輪系角動量與合成標稱角動量的偏差為:

HΣ=CwHw

隨著外擾作用，輪系合成角動量幅值逐步增大直至達到角動量包絡面。當‖HΣ‖≠0時，根據輪系合成角動量實時設定單位向量

S=HΣ/‖HΣ‖

對于由姿態控制律得到的期望控制力矩τc，在偽逆分配律[16]基礎上引入角動量自主調節，得到如下飛輪力矩分配律:

(16)

為對式(16)中所引入的角動量自調整策略特性進行分析，過程中基于如下幾點考慮：

1)空間環境擾動力矩幅值一般均很小，且在短時間內角動量積累量可忽略；

基于上述考慮，在航天器姿態穩定時由式(16)可得：

2)對于由特征值0對應特征向量v0,1，v0,2，v0,3與特征值1對應特征向量v1,1，v1,2，…，v1,N-3所組成的所有特征向量中任意兩者互為正交。

根據上述結論，對于初始時刻t0飛輪角動量Hw(t0)與沿S方向極限角動量的偏差ΔHw(t0)可表示為：

(18)

式(18)右端第一項為常數項，第二項為指數衰減項，并且kNull為ΔHw(t)沿v11，… ，v1,N-3方向指數衰減系數，通過調整該值可改變飛輪角動量調整的動態特性。

CwHw(t)=CwHw(t0)

于是有

即ΔHΣ(t)由角動量包絡上特定點與飛輪當前系統角動量唯一確定?？紤]角動量調節穩態狀態，即在a1,i(t0)=0(i=1,…,N-3)時有

由此可見，上述關于a0,i(t0)(i=1,2,3)為參量的方程組式(19)的解由角動量包絡上特定點與飛輪當前系統角動量唯一確定。

由式(18)形式可知，雖然增大kNull可提高飛輪角動量自調整速度，但當初值ΔHw(t0)過大時可能導致調整所需力矩過大而出現飛輪輸出力矩飽和現象。為克服上述問題以及指數收斂下隨ΔHw調整過程減小使得調整速度變慢問題，除綜合考慮調整速度與調整力矩大小選擇合適kNull值外，可對式(16)角動量調節輸入ΔHw進行單位化，對應飛輪力矩分配律為：

4 仿真與分析

以運行于軌道高度為36 000 km的對地穩定衛星為例，衛星星體轉動慣量為:

符號diag{·}表示對角陣。系統配置4個力矩0.1 N·m和角動量10 N·m·s飛輪，安裝如圖2所示金字塔組合構型。衛星采用整星零動量控制方式，飛輪初始偏置角動量為：

控制參數選取為:

施加于星體的環境外擾力矩為

為對比驗證本文所提出角動量管理策略的效果，首先采用常用的偽逆力矩分配方式進行仿真分析，即設定kNull=0，仿真結果如圖4和圖5所示。結合圖4給出的星體相對軌道系的滾動角φ、俯仰角θ及偏航角Ψ三軸姿態和圖5給出的四飛輪角動量可知，輪系角動量在122 000 s到達角動量包絡面之前飛輪1、飛輪2和飛輪4分別在60 000 s、80 000 s及110 000 s附近不同時間段出現角動量飽和情況，此期間使得姿態角相應地出現一定程度波動。

采用本文提出的角動量管理策略進行仿真，當設定kNull=0.01時在其他相同條件下仿真結果如圖6與圖7所示。結合圖6給出的三軸姿態角和圖7給出的四飛輪角動量可知，輪系在吸收外擾積累角動量過程中通過實時調整飛輪角動量，避免了輪系角動量在到達包絡面之前出現部分飛輪角動量飽和現象且保證控制性能，至到122 000 s兩飛輪角動量同時飽和，即輪系角動量達到包絡面，系統失去三軸姿態控制能力。

在偽逆力矩分配方式下，即kNull=0，增大系統控制參數以分析控制能力提升對系統性能影響。選取控制參數為：

仿真結果如圖8和圖9所示。由圖8可知星體姿態在出現飛輪角動量飽和時姿態波動相比圖4有很大程度改善，但圖9給出的四飛輪角動量與圖5結果基本相同。由上比較可知，雖然通過提高系統反饋控制能力可有效降低飛輪角動量飽和對姿態性能影響，但對系統角動量達到角動量包絡前的飛輪角動量飽和現象無任何改善作用。

通過上述仿真對比可知，采用所提出角動量管理方法可充分利用輪系角動量容納能力，有效避免在輪系角動量積累到達角動量包絡面之前出現部分飛輪飽和情況，相比常規偽逆力矩分配策略可延長高軌衛星噴氣動量卸載間隔。

5 結束語

本文針對受擾輪系航天器角動量管理問題，基于飛輪系統的可視化角動量包絡分析，提出了一種輪系角動量包絡特征參數描述及角動量包絡面處飛輪極限角動量求解方法，在此基礎上給出了一種基于設定角動量目標的零運動力矩分配的角動量管理策略及其特性理論分析，對于既定輪系構型可充分利用輪系角動量空間，可使得飛輪角動量達到角動量包絡面之前避免出現飛輪角動量飽和情況。充分利用輪系角動量的有效性,對所提出力矩分配策略進行了數學仿真驗證。通過與常規偽逆力矩分配策略對比，結果表明該方法可有效延長角動量噴氣卸載間隔，在采用噴氣卸載高軌道衛星角動量管理上具有一定的優勢。

References)

[1] ANDERSEN G C, QUINN D A, BEALS G A, et al. An overview of the Hubble space telescope pointing control system design and operation[C] ∥AIAA Guidance, Navigation, and Control Conference, Hilton Head Island, SC: AIAA, 1992.

[2] 黎康, 喬國棟, 劉新彥. 高分四號衛星控制方案設計特點及在軌實現[J]. 空間控制技術與應用, 2016, 42(6): 1-8.

LI K, QIAO G D, LIU X Y. Highlights of attitude and orbit control system (AOCS) for GaoFen-4 satellite with on-orbit verification [J]. Aerospace Control and Application, 2016, 42(6): 1-8(in Chinese).

[3] 屠善澄. 衛星姿態動力學與控制[M]. 北京:宇航出版社，北京:2001.

[4] 林來興, 白拜爾. 衛星姿態控制反作用輪的最佳安裝結構[J]. 自動化學報，1984, 10(3): 228-238.

LIN L X, BAI B E. The best orientation of reaction wheel in satellite attitude control[J]. Acta Automatic Sinica, 1984, 10(3):228-238(in Chinese).

[5] 顧朋，王大軼，劉成瑞. 零動量衛星輪控系統可重構設計研究[J]. 中國空間科學技術, 2013, 33(1): 7-14.

GU P, WANG D Y, LIU C R. Research on momentum wheel reconfiguration design of zero-momentum satellite [J]. Chinese Space Science and Technology, 2013, 2013, 33(1): 7-14(in Chinese).

[6] 徐帆, 林來興. 角動量包及斜裝飛輪的卸載系統[J]. 航天控制, 1986,4(1): 22-30.

XU F, LIN L X. Angular momentum envelope and unloading system of skewed flywheel system [J]. Aerospace Space Control,1986,4(1): 22-30(in Chinese).

[7] ZULIANA I, RENUGANTH V. A study of reaction wheel configurations for a 3-axis satellite attitude control[J]. Advances in Space Research, 2010, 45(5):750-759.

[8] MARKLEY F L, REYNOLDS R G, LIU F X, et al. Maximum torque and momentum envelopments for reaction wheel arrays[J]. Journal of Guidance, Navigation, and Control Conference, 2010, 33(5):1606-1614.

[9] YOON H, SEO H H, CHOI H T.Optimal uses of reaction wheels in the pyramid configuration using a new minimum infinity-norm solution[J]. Aerospace Science and Technology, 2014, 39: 109-119.

[10] 郭延寧, 馬廣富, 李傳江. 冗余飛輪構型下力矩分配策略設計與分析[J]. 航空學報, 2010, 31(11): 2259-2265.

GUO Y, MA G F, LI C J. Design and analysis of torque allocation strategy for redundant flywheel configuration [J]. Acta Aeronautica ET Astronautic Sinica, 2011, 31(11): 2259-2265(in Chinese).

[11] 黃勇, 胡小工, 董光亮, 等. 動量輪卸載對環月衛星SMART-1衛星軌道的影響和定軌策略[J]. 飛行器測控學報, 2007, 26(5): 62-66.

HUANG Y, HU X G, DONG G L, et al. The effect of wheel off-loading on the orbit of lunar satellite SMART-1[J]. Journal of Spacecraft TT&C Technology, 2007, 26(5): 62-66(in Chinese).

[12] 賴炎連, 賀國平. 最優化方法[M]. 北京:清華大學出版社, 2008:161-216.

[13] 張德豐. Matlab數值分析與應用[M]. 北京:國防工業出版社, 2007:256-260.

[14] 郭雷,王洵,王曉蒲. 有向回路法和網格法:多邊形內外點判別的新算法 [J]. 計算機工程與應用, 2002,38(19):119-122.

GUO L, WANG X, WANG X P.Directed-loop method and grid method: new algorithms for inclusion test of planar polygons[J]. Computer Engineering and Application, 2002, 38(19):119-122(in Chinese).

[15] 韋娟, 寧方立. 誤差四元數及其在航天器姿態控制中的應用[J]. 飛行力學, 2006, 24(2):60-68.

WEI J, NING F L. Error quaternion and its application of attitude control for a spacecraft[J]. Flight Dynamics, 2006, 24(2):60-68(in Chinese).

[16] 章仁為. 衛星軌道姿態動力學與控制[M]. 北京:北京航空航天大學出版社, 1998:273-276.