間接Legendre偽譜法的欠驅(qū)動航天器姿態(tài)運動軌跡跟蹤

易中貴，戈新生

(1. 北京理工大學宇航學院，北京 100081；2. 北京信息科技大學機電工程學院，北京 100192)

0 引 言

近年來，隨著航天科學技術(shù)快速發(fā)展，特別是空間交會對接以及深空探測技術(shù)的發(fā)展使得航天器系統(tǒng)的規(guī)模越來越大，自由度越來越多，構(gòu)造與功能日益復雜，技術(shù)要求逐漸提高，所涉及的領(lǐng)域不斷擴大，需要考慮的因素急劇上升，對運行的速度以及精度也提出了更高的要求。從而增加了姿態(tài)動力學建模、控制仿真分析及運行維護的難度。因此航天器動力學與控制學科的研究受到了越來越多科研工作者的青睞。航天器姿態(tài)控制系統(tǒng)是多個不同功能分系統(tǒng)里面的一個重要分系統(tǒng)。航天器對其姿態(tài)的指向精度有著非常嚴格的要求，如為確保航天器對地通信信號質(zhì)量而使所攜帶的通信天線的方向應長時間精確指向地球方向、航天器展開的太陽帆板需指向太陽方向、空間站等的交會對接過程以及空間機械臂在軌執(zhí)行任務(wù)等等。正常情況下，安裝有足夠或冗余執(zhí)行機構(gòu)(如n個飛輪/控制力矩陀螺，或者n對噴氣推力器(n≥ 3))的三軸穩(wěn)定航天器姿態(tài)控制系統(tǒng)可以在滾轉(zhuǎn)、俯仰和偏航三軸同時輸出控制力矩，完成其姿態(tài)控制和任意定位。但航天器系統(tǒng)的姿態(tài)控制任務(wù)繁多、系統(tǒng)結(jié)構(gòu)復雜，執(zhí)行機構(gòu)長期處于真空、失重、超低溫以及強輻射等實際的惡劣環(huán)境中執(zhí)行在軌任務(wù)容易導致其中的某一個發(fā)生故障或者失效，從而不能輸出完整的三軸控制力矩，故稱此時的航天器系統(tǒng)為欠驅(qū)動航天器。欠驅(qū)動航天器姿態(tài)動力學與控制的研究最早可追溯到20世紀80年代，1984年，Crouch[1]基于微分幾何理論研究了剛體航天器所安裝的執(zhí)行機構(gòu)(噴氣推力器或動量飛輪)數(shù)目少于3個時的可控性問題，并設(shè)計了相應的控制算法，其研究表明，在系統(tǒng)的動量矩不等于零時，欠驅(qū)動航天器并不可控。Krishnan等[2]根據(jù)局部可控性定理研究得到：欠驅(qū)動航天器系統(tǒng)能控的充要條件是其欠驅(qū)動軸不能是其對稱軸。之后國內(nèi)也有不少學者對其進行了更深入的研究和探討，如文獻[3-4]通過動量矩守恒定理先建立了以兩個動量飛輪為執(zhí)行機構(gòu)的欠驅(qū)動航天器的姿態(tài)動力學方程，并在整星零動量的條件下，分別采用Ritz 近似理論以及粒子群優(yōu)化方法，研究了帶兩個反作用動量飛輪的欠驅(qū)動航天器系統(tǒng)的姿態(tài)非完整機動規(guī)劃問題。王芳等[5]研究了航天器上的撓性件以及大氣阻力帶來的正弦干擾力矩情況下的姿態(tài)控制問題中的旋轉(zhuǎn)軸穩(wěn)定問題。鄭敏婕等[6]采用退步控制設(shè)計方法研究了欠驅(qū)動航天器系統(tǒng)的姿態(tài)控制率設(shè)計問題。近年來，張洪華等[7]研究了欠驅(qū)動撓性航天器的姿態(tài)控制問題，提出了“噴氣消旋+飛輪機動”的分段控制方法，并針對此方法的完整過程，采用了擾動系統(tǒng)理論分析閉環(huán)控制系統(tǒng)的大范圍終端有界性。張佳為等[8]針對欠驅(qū)動剛體航天器姿態(tài)控制問題，采用非線性預測控制方法實現(xiàn)了任意飛輪群剪刀對構(gòu)型、飛輪群角動量非飽和條件下，任意系統(tǒng)初始角動量欠驅(qū)動航天器在姿態(tài)可機動集合中的機動控制。宋道喆等[9]通過Lyapunov直接法和Backstepping方法設(shè)計一種非線性不連續(xù)反饋控制律來研究了輪控式零角動量欠驅(qū)動航天器的姿態(tài)最優(yōu)穩(wěn)定控制問題。

對于欠驅(qū)動航天器的姿態(tài)機動控制問題， Legendre偽譜法(Legendre pseudospectral method,LPM)具有較高的精度和較快的運算速度，它通過一系列數(shù)值近似變換，將一個連續(xù)的最優(yōu)控制問題離散為一個非線性規(guī)劃(Non-linear programming，NLP)問題來求解。LPM最初由Elnagar等[10]提出，之后Roa等[11]在此基礎(chǔ)上對其進行分析、擴展，并應用到更廣范的問題上。且Fahroo等[12]的研究提出, 對于一個末端是非完全自由的最優(yōu)控制問題，特別是當末端約束是含初始或終端狀態(tài)的表達式時，Gauss和Radau偽譜法可能不收斂，然而Legendre偽譜法此時卻具有較好的收斂精度和效果。在國內(nèi)，Zhuang等[13]先采用系統(tǒng)的微分平滑特性轉(zhuǎn)化欠驅(qū)動航天器的狀態(tài)和控制輸入變量為一平滑輸出函數(shù)，然后采用LPM離散此平滑函數(shù)為一NLP問題，進而求解出系統(tǒng)的最短時間姿態(tài)機動軌跡。

針對欠驅(qū)動航天器的姿態(tài)運動跟蹤問題，目前出現(xiàn)的研究成果多采用航天器系統(tǒng)的微分平滑特性以及間接LPM對其姿態(tài)運動和跟蹤問題進行研究和探討[14-18]。Aguilar[14]通過證明欠驅(qū)動航天器系統(tǒng)在特殊的轉(zhuǎn)動慣量值時是微分平滑的，然后根據(jù)平滑特性規(guī)劃出系統(tǒng)的開環(huán)參考軌跡，最后再根據(jù)動態(tài)擴展算法設(shè)計跟蹤控制器，從而將線性化系統(tǒng)的姿態(tài)穩(wěn)定到系統(tǒng)的平衡點附近。Tsiotras[15]則先采用系統(tǒng)的微分平滑特性規(guī)劃出系統(tǒng)的可行軌跡，再將其用于后續(xù)的軌跡跟蹤問題中去。文獻[16]基于LPM先規(guī)劃出欠驅(qū)動航天器的最短時間姿態(tài)機動參考軌跡。接著考慮到系統(tǒng)在執(zhí)行任務(wù)過程中不可避免會有初值擾動和不利的干擾激勵等外界因素，進而分別采用間接LPM和系統(tǒng)擴展后的特殊幾何特性設(shè)計了系統(tǒng)的閉環(huán)軌跡跟蹤控制器。文獻[17]則基于改進的Chebyshev偽譜法，采用系統(tǒng)的微分平滑特性以及滾動時域控制方法研究了欠驅(qū)動航天器的姿態(tài)運動軌跡跟蹤控制問題。文獻[18]采用間接LPM研究了高超聲速飛行器再入滑翔的縱向參考軌跡跟蹤制導問題。廖宇新等[19]采用間接Radau偽譜法研究了高超聲速飛行器滑翔段的制導問題，先將標稱軌跡跟蹤問題轉(zhuǎn)化為線性時變系統(tǒng)的狀態(tài)調(diào)節(jié)器問題，然后基于Pontryagin極大值原理將狀態(tài)調(diào)節(jié)器問題轉(zhuǎn)化為線性兩點邊值問題(Two-point boundary value problem,TPBVP)；最后采用間接Radau偽譜法離散所求得的兩點邊值問題為一組線性方程組，避免了在線積分求解Riccati矩陣微分方程等復雜的數(shù)值計算，數(shù)值仿真結(jié)果驗證了所提方法具有良好的魯棒性。

本文研究了欠驅(qū)動航天器存在有較小初始外部擾動時的最短時間姿態(tài)機動優(yōu)化和跟蹤問題。首先采用LPM離線規(guī)劃出系統(tǒng)的最短時間姿態(tài)機動參考軌跡。接著根據(jù)最優(yōu)性原理將跟蹤問題轉(zhuǎn)化為一個兩點邊值問題。鑒于LPM在求解末端約束含有初始或終端狀態(tài)表達式時比Radau偽譜法具有較好的收斂精度[12]。且LPM還具備方陣形式的微分矩陣(離散的線性方程組的個數(shù)等于變量個數(shù))[16]。因此本文基于LPM離散微分方程的思路(并將此方法稱之為間接Legendre偽譜法)，采用Legendre-Gauss-Lobatto (LGL)點離散轉(zhuǎn)化TPBVP為一組線性方程，避免了對傳統(tǒng)Riccati微分方程的積分運算。文末給出了數(shù)值仿真運行結(jié)果。

1 運動模型

根據(jù)歐拉定理可知：空間參考坐標系轉(zhuǎn)動三次即可到達航天器的本體坐標系。本文選取3-2-1的旋轉(zhuǎn)順序。以歐拉角作為姿態(tài)描述參數(shù)的航天器姿態(tài)運動學方程[16]

(1)

假設(shè)航天器系統(tǒng)沒有外力矩作用，且慣量主軸的方向與本體坐標系的方向是同向的，則可寫出航天器的姿態(tài)動力學方程[16]

(2)

(3)

(4)

式中：α=(I1-I2)/I3稱為系統(tǒng)的軸不對稱性系數(shù)。

2 姿態(tài)機動控制

本文先離線規(guī)劃出系統(tǒng)的最短時間姿態(tài)機動軌跡，并將此優(yōu)化軌跡作為后續(xù)姿態(tài)運動軌跡跟蹤的參考軌跡。文中的最短時間控制問題可描述為

minJ=tf=∫tf01dts.tx(t)=f[x(t),u(t),t], t∈[0,tf]x(t0)=x0ψ[x(tf),tf]=0ui∞≤umax(i=1,2)

(5)

其中，ψ[·]表示終端約束條件，umax代表控制輸入力矩的上限值。接著采用Legendre偽譜法離散逼近式(5)所描述的最優(yōu)控制問題，將其轉(zhuǎn)化為一個非線性規(guī)劃問題來求解。鑒于此部分內(nèi)容國內(nèi)外已有較多相關(guān)文獻[10-11]，這里不再贅述。

通過LPM將一個連續(xù)最優(yōu)控制問題轉(zhuǎn)化成的NLP問題可通過SQP算法求解便可得到原時間最優(yōu)控制問題的近似解。文中采用SNOPT軟件包對此NLP問題進行求解。

3 姿態(tài)運動跟蹤

3.1 兩點邊值問題的描述

將欠驅(qū)動航天器系統(tǒng)的動力學方程(式(4))在

已知的參考姿態(tài)運動軌跡上進行泰勒展開，并略去高階項僅保留一次項，從而有線性化方程

(6)

式中：Δx(t)=[Δω1,Δω2,Δω3,Δφ,Δθ,Δψ]T為狀態(tài)變量實際值x與參考值x*的偏差，而Δu(t)=[Δu1,Δu2]T為輸入控制變量u的修正量，Δx0為狀態(tài)初始擾動。時變狀態(tài)矩陣A(t)∈R6×6的表達式為

(7)

時變控制輸入矩陣B(t)∈R6×2的表達式為

(8)

姿態(tài)運動的跟蹤問題可描述為：通過確定控制輸入修正量Δu(t)以及狀態(tài)偏差量Δx(t)，使得如下二次型性能指標最小

ΔuT(t)RΔu(t)]dt

(9)

并滿足線性系統(tǒng)的誤差狀態(tài)方程以及初始條件(式(6))。式中S∈R6×6為半正定的末端加權(quán)矩陣，Q∈R6×6為半正定的狀態(tài)加權(quán)陣，R∈R2×2是正定的控制加權(quán)陣。根據(jù)Pontryagin極小值原理[20]，引入?yún)f(xié)態(tài)變量λ(t)∈R6，并構(gòu)造Hamilton函數(shù)，則最優(yōu)狀態(tài)Δx*(t)和最優(yōu)協(xié)態(tài)λ*(t)滿足正則方程

(10)

邊界條件與最優(yōu)協(xié)態(tài)λ*(t)滿足的橫截條件為

(11)

且根據(jù)Hamilton函數(shù)取極小值時的必要條件有最優(yōu)控制修正量滿足

Δu*(t)=-R-1λ*T(t)B(t)=-R-1BT(t)λ*(t)

(12)

3.2 兩點邊值問題的離散

借鑒LPM離散微分方程的思路。首先通過變換t=(1+τ)tf/2把時間區(qū)間t∈[0,tf]轉(zhuǎn)換到區(qū)間τ∈[-1,1]上。進而由式(10)和式(11)所描述的兩點邊值問題可轉(zhuǎn)化為

(13)

選取K+1個LGL離散配點，將狀態(tài)變量和協(xié)態(tài)變量在這些配點上進行Lagrange插值，并求導有

(14)

將式(14)代入式(13)可得兩點邊值問題的離散表達式為

(15)

Xz=Y

(16)

(17)

(18)

(19)

(20)

至此，整個算法在執(zhí)行過程中不需要任何的積分運算，只需計算線性方程組(16)的解便可求得兩點邊值問題(式(10)～(11))的離散解。此算法既能保證計算精度，又能有效減小計算量和運算時間。

4 數(shù)值仿真

本文先采用LPM離線優(yōu)化出欠驅(qū)動航天器的最短時間姿態(tài)機動軌跡。假設(shè)欠驅(qū)動航天器進行rest-to-rest的姿態(tài)機動控制，其慣量主軸[14]設(shè)為I=diag(66.36,61.80,50.16) kg·m2；初末姿態(tài)約束[14]分別設(shè)置為x0=[0,0,0,-π/2,-π/3,π/2]T，xf=[0,0,0,π/2,π/3,-π/2]T；執(zhí)行機構(gòu)控制輸入上限為umax=1(根據(jù)Pontryagin極小值原理[16,20]，輸入力矩應具有Bang-Bang控制的形式，且在上下限之間切換)；離散點個數(shù)設(shè)置為LGL=37； NLP的求解器為SNOPT[21]，且其求解初值設(shè)為ones(8×38+1,1)。

針對欠驅(qū)動航天器存在初始小擾動進行姿態(tài)跟蹤。先設(shè)置初始偏差值Δx0(見表1)[17]，選擇與參考軌跡相同的LGL配點個數(shù)，代入線性方程組(16)便可得到最優(yōu)協(xié)態(tài)λ*(t)，將其代入方程(12)求得最優(yōu)控制修正量Δu*，進而可得到實際的控制輸入u=u*+Δu*(u*為參考軌跡中規(guī)劃而得的最優(yōu)控制量)。將此控制輸入代入系統(tǒng)運動模型(式(4))進行積分便可求得其實際的姿態(tài)運行軌跡x。狀態(tài)偏差可由Δx=x-x*得到。加權(quán)矩陣S，Q和R可根據(jù)Bryson’s Rule確定[22]

(21)

式中：下標max表示各變量所允許的最大偏差，本文取狀態(tài)變量的最大偏差為2Δx0，而控制變量的最大偏差取為[0.0055, 0.0055]。并設(shè)Q11=1，則矩陣Q和R中的其余元素便可由式(21)確定。

具體仿真結(jié)果見圖1～圖5，圖中同時給出了姿態(tài)機動參考軌跡以及兩種極限初始小擾動情況下的閉環(huán)和開環(huán)跟蹤控制曲線。

表1 初始擾動Table 1 The initial disturbance

注:表1中Δω的單位為rad/s，三個歐拉角的單位為度。

圖1 參考及跟蹤控制輸入Fig.1 The reference and tracking control torques

圖2 角速度的參考、開環(huán)以及閉環(huán)軌跡Fig.2 The trajectories for angular velocities

圖3 姿態(tài)角的參考、開環(huán)以及閉環(huán)軌跡Fig.3 The trajectories for Euler angles

圖4 角速度跟蹤誤差曲線Fig.4 The tracking errors for angular velocities

圖5 姿態(tài)角跟蹤誤差曲線Fig.5 The tracking errors for Euler angles

從圖1中的控制輸入滿足開關(guān)控制(Bang-Bang控制)的形式，且在預設(shè)的最大和最小控制值之間切換。從圖2～3中的參考軌跡可以看出，三個歐拉姿態(tài)角和姿態(tài)角速度均可精確規(guī)劃到預定的末端值，優(yōu)化曲線均平滑，優(yōu)化目標的最短時間為5.0166 s。

圖2～3中的點劃線和虛線分別代表上下極限初始擾動情況下的閉環(huán)跟蹤控制曲線，而細實線和點線則分別為帶有上下極限初始小擾動的開環(huán)跟蹤控制曲線。從圖2～3可以看出，閉環(huán)跟蹤效果較好，對初始外部擾動具有一定的抑制能力，而開環(huán)則不能跟蹤到目標點，這可從圖4～5中的跟蹤誤差曲線得到體現(xiàn)。

5 結(jié) 論

本文研究了帶有兩組噴氣推力器的非對稱欠驅(qū)動航天器的姿態(tài)優(yōu)化控制及其姿態(tài)運動跟蹤問題。先采用LPM規(guī)劃出欠驅(qū)動航天器的最短時間姿態(tài)機動參考軌跡。接著采用LGL點離散時變偏差系統(tǒng)的TPBVP，將其轉(zhuǎn)化為一組線性方程組來求解，避免了對Riccati方程的積分運算。最后通過數(shù)值仿真校驗了本文所提出算法的有效性和合理性。

參 考 文 獻

[1] Crouch P E. Spacecraft attitude control and stabilization: applic-ation of geometric control theory to rigid body models[J]. IEEE Transactions on Automatic Control, 1984, 29(4): 87-95.

[2] Krishnan H, Mc clamroch N H, Reyhanoglu M. Attitude stabilization of a rigid spacecraft using two momentum wheel actuators[J]. Journal of Guidance Control & Dynamics, 1995, 18(2): 256-263.

[3] 戈新生, 陳立群, 劉延柱. 帶有兩個動量飛輪剛體航天器的姿態(tài)非完整運動規(guī)劃問題[J]. 控制理論與應用, 2004, 21(5): 781-784. [Ge Xin-sheng, Chen Li-qun, Liu Yan-zhu. Nonholonomic motion planning for the attitude of rigid spacecraft with two momentum wheel actuators[J]. Control Theory & Applications, 2004, 21(5): 781-784.]

[4] 戈新生, 孫鵬偉. 欠驅(qū)動剛性航天器姿態(tài)的非完整運動規(guī)劃粒子群算法[J]. 宇航學報, 2006, 27(6): 1233-1237. [Ge Xin-sheng, Sun Peng-wei. Nonholonomic motion planning for the attitude of underactuated spacecraft using particle swarm optimization[J]. Journal of Astronautics, 2006, 27(6): 1233-1237.]

[5] 王芳, 張洪華. 欠驅(qū)動剛性航天器旋轉(zhuǎn)軸穩(wěn)定研究[J]. 宇航學報, 2007, 28(5):1133-1137. [Wang Fang, Zhang Hong-hua. Spin-axis stabilization of underactuated rigid spacecraft[J]. Journal of Astronautics, 2007, 28(5):1133-1137.]

[6] 鄭敏捷, 徐世杰. 欠驅(qū)動航天器姿態(tài)控制系統(tǒng)的退步控制設(shè)計方法[J]. 宇航學報, 2006, 27(5):63-66. [Zheng Min-jie, Xu Shi-jie. Backstepping control for attitude control system of an underactuated spacecraft[J]. Journal of Astronautics, 2006, 27(5):63-66.]

[7] 張洪華, 王芳, 胡錦昌, 等. 欠驅(qū)動撓性航天器的全姿態(tài)控制[J]. 宇航學報, 2015, 36(4):419-429. [Zhang Hong-hua, Wang Fang, Hu Jin-chang, et al. All-attitude control for underactuated flexible spacecraft[J]. Journal of Astronautics, 2015, 36(4):419-429.]

[8] 張佳為, 許諾, 伍少雄. 欠驅(qū)動航天器飛輪控制方法[J]. 宇航學報, 2016, 37(5):552-561. [Zhang Jia-wei， Xu Nuo， Wu Shao-xiong. Control method of underactuated spacecraft using flywheels[J]. Journal of Astronautics, 2016, 37(5):552-561.]

[9] 宋道喆, 耿云海, 易濤. 零角動量欠驅(qū)動航天器逆最優(yōu)穩(wěn)定控制律設(shè)計[J]. 宇航學報, 2016, 37(6):729-736. [Song Dao-zhe， Geng Yun-hai， Yi Tao. Inverse optimal stabilization of an underactuated spacecraft in a zero angular momentum mode[J]. Journal of Astronautics, 2016, 37(6):729-736.]

[10] Elnagar G N, Kazemi M A. The pseudospectral Legendre method for discretizing optimal control problems[J]. IEEE Transactions on Automatic Control, 1995, 40(10):1793-1796.

[11] Rao A V, Clarke K A. Performance optimization of a maneuvering reentry vehicle using a Legendre pseudospectral method[J]. AIAA Journal, 2002.

[12] Fahroo F, Ross I M. Advances in pseudospectral methods for optimal control[C]. AIAA Guidance, Navigation and Control Conf and Exhibit, Honolulu,USA, 2008.

[13] Zhuang Y F, Guang-Fu M A, Chuan-Jiang L I, et al. Time-optimal trajectory planning for underactuated rigid spacecraft using differential flatness[J]. Journal of Astronautics, 2011, 32(8):1753-1761.

[14] Aguilar C O. Attitude control of a differentially flat underactuated rigid spacecraft [D]. Edmonton: University of Alberta, 2005.

[15] Tsiotras P. Feasible trajectory generation for underactuated spacecraft using differential flatness[J]. Advances in the Astronautical Sciences, 2008, 103:16-18.

[16] 莊宇飛. 帶有非完整約束的欠驅(qū)動航天器控制方法研究[D]. 哈爾濱: 哈爾濱工業(yè)大學, 2012. [Zhuang Yu-fei. Control of underactuated spacecraft with nonholonomic constraints[D]. Harbin: Harbin Institute of Technology, 2012.]

[17] Cai W W, Yang L P, Zhu Y W. Optimal reorientation of asymmetric underactuated spacecraft using differential flatness and receding horizon control[J]. Advances in Space Research, 2015, 55(1):343-353.

[18] 李強. 高超聲速滑翔飛行器再入制導控制技術(shù)研究[D]. 北京: 北京理工大學, 2015. [Li Qiang. Study on reentry guidance and control method for hypersonic glide vehicle[D]. Beijing: Beijing Institute of Technology, 2012.]

[19] 廖宇新, 李惠峰, 包為民. 基于間接Radau偽譜法的滑翔段軌跡跟蹤制導律[J]. 宇航學報, 2015, 36(12):1398-1405. [Liao Yu-xin, Li Hui-feng, Bao Wei-min. Gliding trajectory tracking guidance law based on indirect Radau pseud-ospectral method[J]. Journal of Astronautics, 2015, 36(12):1398-1405.]

[20] 李傳江，馬廣富. 最優(yōu)控制[M]. 北京: 科學出版社, 2011.

[21] Gill P E, Murray W, Saunders M A. SNOPT: an SQP algorithm for large-scale constrained optimization[J]. Siam Journal on Optimization, 2002, 12(1):99-131.

[22] Bryson D E, Ho Y C. Applied optimal control[M]. Washington, D C: Hemisphere Publishing Corporation, 1975.