蔡偉偉 朱彥偉 曾璞

(1.國防科技大學(xué)空天科學(xué)學(xué)院,長沙 410073) (2.國防科技大學(xué)軍事基礎(chǔ)教育學(xué)院,長沙 410073)

引言

空間飛行器在軌運行期間,需要通過大量的姿態(tài)機動在不同姿態(tài)指向之間進行切換,達成多樣化的任務(wù)目的；然而,姿態(tài)執(zhí)行機構(gòu)性能、星上傳感器空間指向等復(fù)雜約束的存在,使得姿態(tài)機動路徑面臨諸多限制,需要在多約束條件下開展姿態(tài)機動路徑規(guī)劃研究,為安全、穩(wěn)定地工作提供支撐.針對空間飛行器姿態(tài)機動規(guī)劃問題,國內(nèi)外學(xué)者開展了許多有意義的研究.武長青從能量最優(yōu)角度出發(fā),將航天器姿態(tài)約束機動問題歸納為非凸二次約束二次規(guī)劃問題,運用基于評價函數(shù)的迭代規(guī)劃算法獲得全局姿態(tài)優(yōu)化路徑[1].趙乾考慮力矩陀螺奇異與飽和特性,根據(jù)工程實際偏好需求,研究了空間站姿態(tài)機動路徑規(guī)劃的多目標(biāo)優(yōu)化問題[2].豐志偉針對單框架控制力矩陀螺驅(qū)動下的敏捷衛(wèi)星姿態(tài)機動問題,考慮奇異狀態(tài)規(guī)避,利用偽譜法對時間固定、能量最優(yōu)的姿態(tài)機動路徑進行了優(yōu)化[3].

近年來,一種基于微分平坦理論的軌跡規(guī)劃方法因其顯著的計算效率而廣泛應(yīng)用于四旋翼飛行器[4]、動態(tài)滑翔無人機[5]等領(lǐng)域.該方法利用系統(tǒng)的微分平坦屬性,將初始軌跡規(guī)劃問題映射到平坦輸出空間,有效消除微分約束,降低設(shè)計維度,從而達到提高計算效率的目的.在空間飛行器姿態(tài)機動方面,莊宇飛針對欠驅(qū)動航天器的姿態(tài)機動規(guī)劃問題,引入虛擬控制輸入的概念擴展了系統(tǒng)微分平坦屬性,并利用樣條函數(shù)、三角多項式函數(shù)等參數(shù)化平坦輸出,實現(xiàn)姿態(tài)機動可行路徑的快速規(guī)劃[6].然而,函數(shù)逼近論表明上述參數(shù)化方法在近似精度和計算效率兩方面均不是最佳的選擇.因此,對于最優(yōu)軌跡規(guī)劃問題需尋求其他平坦輸出參數(shù)化方法.由于偽譜法能夠以較少的離散點對待定函數(shù)進行高精度近似,Ross最早將偽譜法應(yīng)用于平坦輸出參數(shù)化,并直接由偽譜微分矩陣計算平坦輸出的各階導(dǎo)數(shù),改善了微分平坦方法的性能[7].然而,標(biāo)準(zhǔn)偽譜法的微分矩陣存在病態(tài)特性,嚴重影響平坦輸出高階導(dǎo)數(shù)在規(guī)劃時域端點處的精度.因此,為提高機動軌跡規(guī)劃的性能,需進一步改善偽譜微分矩陣的病態(tài)特性.

本文研究了復(fù)雜約束條件下空間飛行器姿態(tài)機動的規(guī)劃問題,建立了姿態(tài)機動軌跡軌跡規(guī)劃的模型,分析了姿態(tài)運動方程的微分平坦屬性,分別給出了基于微分平坦的規(guī)劃策略和映射Chebyshev偽譜法參數(shù)化平坦輸出的方法,并開展數(shù)值仿真分析驗證所研究方法的性能.

1 姿態(tài)機動規(guī)劃問題描述

1.1 姿態(tài)運動模型

本文以非對稱剛體空間飛行器為研究對象,其沿各慣量主軸方向的姿態(tài)動力學(xué)方程為:

(1)

采用歐拉角描述空間飛行器姿態(tài),按3-2-1轉(zhuǎn)動順序,則姿態(tài)運動學(xué)方程為:

(2)

式中φ,θ和ψ分別為滾轉(zhuǎn)、俯仰和偏航角.

1.2 約束條件

在初始和終端期望姿態(tài)之外,空間飛行器姿態(tài)還面臨控制輸入受限、傳感器指向受限等復(fù)雜約束.實際工程中,姿態(tài)執(zhí)行機構(gòu)提供的力矩幅值有限,形成控制輸入飽和約束:

ui(t)≤umax,t∈[t0,tf]

(3)

式中i=1,2,3;t0和tf分別表示初始和終端時刻.

在姿態(tài)機動過程中,紅外敏感望遠鏡等傳感器不能朝向太陽,以免損害敏感元件,這類約束稱為空間指向約束:

(4)

式中rB表示傳感器在空間飛行器本體系下的方向矢量,rI表示太陽在慣性系下的方向矢量,CBI為慣性系相對于本體系之間的姿態(tài)矩陣,β為rB和rI間的最大夾角.

1.3 目標(biāo)函數(shù)

在達成期望姿態(tài)的同時,姿態(tài)機動過程還可對機動時間、能量等進行優(yōu)化,本文以能量最優(yōu)為控制目標(biāo):

(5)

2 姿態(tài)機動規(guī)劃方法

2.1 微分平坦簡介

對于一般非線性系統(tǒng):

(6)

其中x∈n為狀態(tài)變量,u∈m為控制輸入,f:n×m→n為光滑函數(shù).若存在一組相互微分獨立的變量ζ=[ζ1,ζ2…ζm]∈m滿足:

(7)

且系統(tǒng)狀態(tài)和輸入能夠表示為:

(8)

(9)

則稱非線性系統(tǒng)(6)是微分平坦的,其中變量ζ為平坦輸出[8].

根據(jù)微分平坦系統(tǒng)定義,平坦輸出與系統(tǒng)狀態(tài)和控制輸入之間存在可逆映射關(guān)系,系統(tǒng)運動行為能夠完全由平坦輸出決定.基于此,可將狀態(tài)空間中的軌跡規(guī)劃問題映射到平坦輸出空間中,解出最優(yōu)平坦輸出函數(shù)后,重新映射回狀態(tài)空間即可生成期望的狀態(tài)和控制變量軌跡.

2.2 姿態(tài)運動模型微分平坦屬性

對第1.1節(jié)中的姿態(tài)運動模型,以歐拉角φ,θ和ψ為候選平坦輸出,分別記為ζ1=φ,ζ2=θ,ζ3=ψ,則其余狀態(tài)變量可表示為:

(10)

將式(10)對時間t求導(dǎo)得:

(11)

將式(10)和(11)代入式(1)中,則控制輸入變量可表示為平坦輸出及其各階導(dǎo)數(shù)的函數(shù):

(12)

至此可以判斷姿態(tài)運動模型具有微分平坦屬性,將式(10)和(12)代入約束條件式(3)～(4)和目標(biāo)函數(shù)式(5)中,則姿態(tài)機動規(guī)劃模型可以映射到平坦輸出空間中:確定平坦輸出ζ=[ζ1,ζ2,ζ3]T= [φ,θ,ψ]T,使得如下目標(biāo)函數(shù)最小:

(13)

滿足邊界條件、控制輸入飽和與空間指向約束:

(14)

(15)

(16)

其中式(14)和(16)分別為平坦輸出表示的邊界條件與空間指向約束.

由于平坦輸出是微分平坦系統(tǒng)的最小描述,其維數(shù)與非線性系統(tǒng)控制輸入維數(shù)相等,因此對于狀態(tài)空間維數(shù)較高的系統(tǒng),基于微分平坦理論將軌跡規(guī)劃問題映射到平輸出空間中,能夠有效減少系統(tǒng)分析設(shè)計的維度,降低問題求解難度.此外,通過微分平坦變換,系統(tǒng)微分方程約束得以完全消除,轉(zhuǎn)換所得幾何規(guī)劃問題更易求解,計算效率更高.

2.3 映射Chebyshev偽譜法

本節(jié)給出參數(shù)化平坦輸出的映射Chebyshev偽譜法,該方法以映射Chebyshev-Gauss-Lobatto(CGL)節(jié)點為插值節(jié)點,運用重心有理Lagrange插值技術(shù)近似平坦輸出,由微分矩陣直接獲取平坦輸出在配點處的導(dǎo)數(shù),并通過Clenshaw-Curtis積分替代目標(biāo)函數(shù)中的積分項,最終將平坦輸出規(guī)劃問題進一步轉(zhuǎn)換為非線性規(guī)劃問題求解.相對于高階多項式、樣條插值等方式,該方法能夠以較少的節(jié)點獲得平坦輸出高精度近似；相對于標(biāo)準(zhǔn)Chebyshev偽譜法,該方法能有效改善微分矩陣的病態(tài)特性,有利于提高平坦輸出高階導(dǎo)數(shù)的近似精度[9].

(1)時域變換

映射Chebyshev偽譜法的插值節(jié)點分布在區(qū)間[-1,1]上,因此通過式(17)所示變換將姿態(tài)機動的時間區(qū)間映射到閉區(qū)間[-1,1]上:

λ=2(t-t0)/(tf-t0)-1,t∈[t0,tf]

(17)

(2)平坦輸出近似

映射CGL點λk(k=0,…,N)定義為:

λk=g(τk,α)=asin(ατk)/asinα

(18)

式中0≤α<1為Kosloff-Tal-Eaer共形映射參數(shù),τk為第Ⅱ類標(biāo)準(zhǔn)CGL節(jié)點:

τk=cos(πk/N)

(19)

共形映射的一對一以及充分光滑特性有效保持了Chebyshev偽譜法的收斂特性,且共形映射引入的舍入誤差ε與參數(shù)α、N間存在下列關(guān)系:

α=sech(|lnε|/N)

(20)

當(dāng)舍入誤差ε與機器誤差一致時,可有效降低共形映射引入的誤差.

值得注意的是,映射CGL點相對于標(biāo)準(zhǔn)CGL點更趨于均勻分布；若直接采用傳統(tǒng)Lagrange插值技術(shù)近似平坦輸出,將可能產(chǎn)生龍格現(xiàn)象,因此,此處采用重心有理Lagrange插值技術(shù).記平坦輸出ζi在映射CGL點處的取值為ζi(λk)(i=1…m;k=0,1…N),則平坦輸出近似表達式為:

(21)

(22)

(3)平坦輸出導(dǎo)數(shù)近似

平坦輸出在映射CGL點處的導(dǎo)數(shù)也可以寫成微分矩陣乘積的形式:

(23)

若i≠j,

(24)

若i=j,

(25)

(4)目標(biāo)函數(shù)近似

映射Chebyshev偽譜法采用Clenshaw-Curtis積分公式轉(zhuǎn)換式(13)所示的積分型性能指標(biāo):

(26)

若N為偶數(shù),則:

(27)

式中s=1,…,N/2.

若N為奇數(shù),則:

(28)

式中s=1,…,(N-1)/2,求和符號上方的兩撇表示首末兩項表達式應(yīng)除以2.

綜上,采用映射Chebyshev偽譜法參數(shù)化平坦輸出,微分平坦變換所得最優(yōu)問題式(13)～(16)可以轉(zhuǎn)換為如下非線性規(guī)劃問題:確定平坦輸出在映射CGL點處的取值ζ(λk)(k=0…N)以及終端時刻tf,使得式(26)所示性能指標(biāo)最小,同時滿足邊界條件、控制輸入飽和與空間指向約束:

B(ζ(-1),t0;ζ(1),tf)=0

(29)

Γui(ζ(λk),t)≤umax,(i=1,2,3)

(30)

C(ζ(λk),t)≤Cmax

(31)

通過序列二次規(guī)劃等數(shù)值方法求得平坦輸出最優(yōu)解后,代入狀態(tài)和控制變量關(guān)于平坦輸出的表達式,即式(10)和(12)中,即可獲得姿態(tài)機動規(guī)劃問題的最優(yōu)解.

3 仿真分析

通過微分平坦理論對問題進行處理,并采用映射Chebyshev偽譜法參數(shù)化平坦輸出,其中Chebyshev多項式階次N和舍入誤差ε分別取為32和1.0×10-6,則共形映射參數(shù)α為0.9135,共有決策變量99個.對最終轉(zhuǎn)換所得非線性規(guī)劃問題,調(diào)用MATLAB?軟件中的fmincon函數(shù)進行求解,其參數(shù)TolFun、TolX和TolCon均設(shè)置為1.0×10-6.對于決策變量的初始猜測值,即平坦輸出在映射CGL點處取值,通過樣條插值邊界條件獲得.fmincon函數(shù)經(jīng)246次迭代計算,規(guī)劃所得能量最省機動路徑對應(yīng)的歐拉角、角速度和控制輸入分別如圖1～3所示,圖中圓圈表示相應(yīng)變量在映射CGL點處的取值,實線為經(jīng)重心有理Lagrange插值獲得的相應(yīng)變量隨時間變化關(guān)系.式(5)所示目標(biāo)函數(shù)取值為0.2607.

圖1 歐拉角-時間曲線Fig.1 Time histories of Euler angles

為驗證微分平坦方法性能,基于開源優(yōu)化軟件GPOPS,運用Radau偽譜法對相同配置的問題進行求解,其中配點數(shù)同樣為32,計算網(wǎng)格固定[10].值得說明的是,GPOPS調(diào)用商業(yè)軟件包SNOPT求解Radau偽譜法轉(zhuǎn)換所得非線性規(guī)劃問題.兩種方法的性能對比如表1所示,顯然微分平坦方法有效減少了決策變量數(shù)量,且迭代優(yōu)化的次數(shù)大幅降低.運用GPOPS規(guī)劃所得軌跡與微分平坦法規(guī)劃結(jié)果的差異水平如圖4～6所示,兩者差值的數(shù)量級達到1.0×10-4以上.由此可見,基于微分平坦的姿態(tài)機動規(guī)劃方法有效降低了規(guī)劃空間的維度,有利于規(guī)劃效率的提高.

圖2 角速度-時間曲線Fig.2 Time histories of angle velocity

圖3 控制輸入-時間曲線Fig.3 Time histories of control inputs

圖4 歐拉角差異Fig.4 Errors of the Euler angles

MethodNumber of decision variablesNumber of iterationsObjective functionRadau2944700.2607Proposed method992460.2607

圖5 角速度差異Fig.5 Errors of the angle velocity

圖6 控制輸入差異Fig.6 Errors of the control inputs

4 小結(jié)

本文研究了控制輸入飽和、空間指向受限等多約束條件下的空間飛行器姿態(tài)機動規(guī)劃問題,提出基于微分平坦理論的規(guī)劃方法,將問題映射到平坦輸出空間,在消除微分約束的同時降低了規(guī)劃空間的維度；給出了平坦輸出參數(shù)化的映射Chebyshev偽譜法,最終轉(zhuǎn)換為非線性規(guī)劃問題求解.數(shù)值仿真驗證了本文方法的計算性能,對工程應(yīng)用具有一定參考價值.值得說明的是,文中采用歐拉角描述空間飛行器姿態(tài),其在大角度姿態(tài)機動情況下存在奇異問題,因此本文方法解決以四元數(shù)、修正羅德里格斯參數(shù)等描述的姿態(tài)機動規(guī)劃問題是需要進一步研究的內(nèi)容.

1武長青,徐瑞. 非凸二次約束下航天器姿態(tài)機動路徑迭代規(guī)劃方法. 宇航學(xué)報, 2016,37(6):671～678 (Wu C Q, Xu R. Spacecraft attitude maneuver path iterative planning method under nonconvex quadratic constraints.JournalofAstronautics, 2016,37(6):671～678 (in Chinese))

2趙乾. 考慮力矩陀螺奇異與飽和的空間站姿態(tài)機動控制策略與路徑規(guī)劃[博士學(xué)位論文]. 長沙:國防科技大學(xué), 2015 (Zhao Q. Study on space station attitude maneuver control strategy and path planning considering saturation and singularity of control moment gyroscopes[Ph.D Thesis]. Changsha: National University of Defense Technology, 2015(in Chinese))

3豐志偉,張永合. 基于路徑規(guī)劃的敏捷衛(wèi)星姿態(tài)機動反饋控制方法. 國防科技大學(xué)學(xué)報, 2013,35(4):1～6 (Feng Z W, Zhang Y H. Feedback control method for attitude maneuver of agile satellite based on trajectory optimization.JournalofNationalUniversityofDefenseTechnology, 2013,35(4):1～6(in Chinese))

4Chamseddine A, Zhang Y, Rabbath C A, et al. Flatness-based trajectory planning/replanning for a quadrotor unmanned aerial vehicle.IEEETransactionsonAerospaceandElectronicSystems, 2012,48(4):2832～2848

5Deittert M, Richards A, Toomer C A, et al. Engineless unmanned aerial vehicle propulsion by dynamic soaring.JournalofGuidance,ControlandDynamics, 2009, 32(5): 1446～1457

6莊宇飛. 帶有非完整約束的欠驅(qū)動航天器控制方法研究[博士學(xué)位論文]. 哈爾濱:哈爾濱工業(yè)大學(xué), 2012 (Zhuang Y F. Control of underactuated spacecraft with nonholonomic constraints[Ph.D Thesis]. Harbin: Harbin Institute of Technology, 2012)

7Ross I M, Fahroo F. Pseudospectralmethods for optimal motion planning of differentially flat systems.IEEETransactionsonAutomaticControl, 2004,49(8):1410～1413

8Fliess M. Flatness and Defect of Nonlinear Systems: Intorductory Theory and Examples.InternationalJournalofControl, 1995,61(6):1327～1361

9Cai W W, Zhu Y W, Yang L P. Optimal guidance for hypersonic reentry using inversion and receding horizon control.LetControlTheory&Applications, 2015,9(9):1347～1355

10 Rao A V, Benson D A, Darby C, et al. Gpops, a Matlab Software for Solving Multiple-Phase Optimal Control Problems Using the Gauss Pseudospectral Method.ACMTransactionsonMathematicalSoftware, 2010,37(2):1～39