馮逸駿, 陳萬春, 楊良
(北京航空航天大學 宇航學院, 北京 100083)
現代空間任務通常要求飛行器具有良好的姿態機動能力[1]。空間飛行器的姿態機動是指飛行器在外太空作大角度的姿態調整[2]。相比于人造衛星、空間站等航天器,空間飛行器(如亞軌道飛行器、動能攔截器、軌道再入攻擊器等,以下簡稱飛行器)的姿態機動存在以下特點:①姿態機動角度大、時間短;②一般為主動控制,多采用反作用控制系統(Reaction Control Systems, RCS),甚至是定推力的RCS進行控制;③關心能量消耗,在滿足精度條件下要求機動耗能盡可能小;④由于機動時間短,通常忽略如重力梯度、太陽光壓等擾動。以上特點為飛行器姿態機動控制系統的設計提出了新的挑戰。
飛行器在外太空姿態機動的過程中通常可看作為剛體。描述六自由度剛體姿態的系統通常是具有強耦合特性的多輸入多輸出(MIMO)非線性系統,因此飛行器能量最優姿態機動控制問題實質上是有限時域內非線性系統最優控制問題的一個特例。 針對此問題,Bharadwaj等提出了逆最優控制(inverse optimal control)的方法[3],通過構造李雅普諾夫函數和求解HJB(Hamilton-Jacobi-Bellman)方程獲得反饋控制量,實現滿足性能指標最優的動力學系統控制,并被廣泛用于飛行器能量最優姿態機動和控制當中,如文獻[4-7]。然而,逆最優控制方法滿足的是終端精度與能量消耗的加權和在無窮時域上的最優,并非有限時間內滿足終端條件的能量最優控制。同時,逆最優控制姿態機動的終端精度、機動時間和消耗能量取決于若干控制參數的調節,在實際使用中有所不便。
為此,本文提出了一種基于模型預測控制(Model Predictive Control,MPC)和線性偽譜的能量最優姿態機動控制方法。
模型預測控制是一種基于滾動優化的在線控制策略,具有對模型要求低、抗干擾性好、魯棒性強等優點,能夠在優化性能指標的同時處理各種約束條件,得到了工程技術人員和理論研究者的重視,并被應用到飛行器姿態控制問題的研究當中[1,8-9]。近年來,快速MPC逐漸引起關注。模型預測控制通常涉及到當前狀態線性化后的局部最優控制問題,一般地,求解該最優控制問題可轉化為求解一個兩點邊值問題(Two-Point Boundary Value Problem,TPBVP),專家學者為提高控制指令求解速率和控制精度做了大量工作[10],針對具有強終端約束以及二次性能指標的控制問題,提出了模型預測靜態規劃(Model Predictive Static Programming,MPSP)非線性最優控制方法[11],并將其成功運用在飛行器再入制導領域[12]。盡管MPSP方法能獲得全局最優解,但其在計算過程中需要采用大量離散節點以保證積分精度,導致計算效率較低。針對該問題,Yang等[13]提出了線性高斯偽譜模型預測控制(Linear Gauss Pseudo-spectral Model Predictive Control,LGPMPC)方法,該方法綜合了非線性近似模型預測控制、線性二次最優控制以及高斯偽譜法,采用較MPSP更少的離散節點達到更高的計算效率及精度,被應用于大氣層外制導[14]。針對多段問題,Yang等[15]又在LGPMPC的基礎上提出了多段線性偽譜模型預測控制(Multi-segment Linear Pseudo-spectral Model Predictive Control,MLPMPC)方法,并被成功應用于飛行器再入段制導。
本文在LGPMPC和MLPMPC方法的基礎上,針對飛行器姿態機動問題,利用KKT (Karush-Kuhn-Tucher)條件和線性偽譜推導了能量最優控制指令修正量的解析表達式,從而得到基于線性偽譜模型預測控制的能量最優姿態機動控制方法。該方法能夠同時滿足多個終端約束,并保證修正后的控制指令仍然滿足能量最優。基于線性偽譜模型預測控制的能量最優姿態機動控制方法的基本策略為:根據飛行任務離線優化出初始姿態機動軌跡,實際飛行過程中在線修正標稱控制以保證姿態機動終端精度,同時迭代更新姿態機動標稱軌跡。數值仿真結果表明,該方法能夠在限定時間內實現能量最優大角度姿態機動,在近似精度的前提下,相比于傳統線性二次型調節器(Linear Quadratic Regulator,LQR)跟蹤節省約10%的能量消耗。
本文采用修正羅德里格斯參數(Modified Rodrigues Parameters,MRPs)描述飛行器姿態。考慮一般的剛體飛行器,其姿態運動學和動力學方程可描述為

(1)

(2)
G(σ)∈R3×3為飛行器姿態運動學矩陣,定義為
(3)
式中:I3×3為3×3單位矩陣。飛行器姿態機動過程中,控制力矩存在如式(4)幅值約束:
|ui|≤Uimaxi=1,2,3
(4)
式中:Uimax>0為每個姿態控制通道上的最大力矩幅值。
設狀態量x=[σ1σ2σ3ω1ω2ω3]T,則飛行器姿態運動系統可采用狀態空間的形式描述為
(5)



(6)
式中:t0和tf分別為飛行器姿態機動的起始時刻和終端時刻;φ為末值型項指標;φ為積分型項指標。此處性能指標Φ只考慮能量最優,故只與控制量相關,一般地,取
(7)
線性偽譜模型預測控制方法的主要思想為:將非線性動力學方程在標稱狀態量附近進行擬線性化,建立一個以偏差為自變量的線性微分方程,通過拉格朗日插值多項式對該線性微分方程的狀態量進行逼近,將微分動力學約束通過正交配點轉化為一組多變量代數約束,最終將非線性方程的積分問題轉化為一個連續求解線性代數方程組的問題,從而得到非線性動力學過程的終端狀態與狀態偏差、控制偏差的等式關系,以實現偏差修正[13-14]。
基于線性偽譜模型預測控制的能量最優姿態機動控制策略需要基于一條標稱機動軌跡,因此須先求解最優控制問題,規劃得到標稱機動軌跡。本文采用高斯偽譜法將能量最優姿態機動問題離散為非線性規劃問題,采用SNOPT工具包求解該非線性規劃問題。標稱機動軌跡的具體求解及優化過程,并非本文研究重點,此處不再贅述。
線性偽譜模型預測能量最優姿態機動控制方法可以分為終端狀態解析預測和控制指令修正兩部分。
一般地,考慮具有終端約束的姿態機動非線性動力學方程如下:

(8)


(9)
(10)


(11)
式中:A(t)∈R6×6為狀態誤差傳播矩陣,B(t)∈R6×3為控制誤差傳播矩陣。對式(1)~式(3)描述的飛行器姿態動力學系統而言,狀態誤差傳播矩陣和控制誤差傳播矩陣的具體形式為
(12)
其中:A11∈R3×3,A12∈R3×3,A22∈R3×3,具體表達式分別為
(13)
(14)
(15)
式中:K、S1、S2和S3為中間變量,其表達式分別為
(16)
(17)
(18)
(19)
對狀態量和控制量進行偽譜離散時,通常有3種高斯正交節點可供選擇:LG(Legendre Gauss)節點、LGR(Legendre Gauss Radau)節點、LGL(Legendre Gauss Lobatto)節點。其中,LG節點不含終端點,LGR節點含有一個終端點(通常為左端歸一化時間τ=-1),LGL節點含有2個終端點,如圖1所示。考慮到工程實際中的控制連續性,此處選取LGR節點對狀態量和控制量進行離散。
首先,將實際機動時間t∈[t0,tf]映射到歸一化時間τ∈[-1,1]上:
(20)
則歸一化后的誤差傳播動力學方程為
(21)
式中:p為歸一化轉換變量,其表達式為
(22)

(23)
同樣地,本文得到控制量的插值擬合:
(24)
拉格朗日插值多項式滿足以下性質:
(25)
故有
(26)
通過對狀態量求導可得
(27)
式中:微分逼近矩陣D∈RN×(N+1)是通過對拉格朗日插值多項式的各個元素分別求導獲得的,其具體表達式為

圖1 高斯正交節點示意圖Fig.1 Schematic diagram of Gauss quadrature points

(28)

(29)
式中:k=1,2,…,N。

(30)
重新組合微分逼近矩陣D,可得到預測狀態偏差序列和控制修正序列表示的關系式如下:
(31)
其中:
(32)
(33)
(34)
矩陣A*和B*的表達式分別為
(35)
(36)
則其他LGR節點上的各個狀態量可表示為
(37)

(38)


(39)
其中:Wx為高斯型積分公式的權函數矩陣,可表示為
(40)
(41)
(42)
式中:Kx∈R6×6,Ku∈R6×3N,具體表達式分別為
Kx=w0pA(τ0)-
pA′]-1D1
(43)
pA′]-1B*
(44)
至此,本文獲得了預測終端狀態偏差關于初始狀態偏差和控制修正的解析等式關系。
在求得預測終端狀態偏差關于初始狀態偏差和控制修正的解析等式關系后,可以根據最優控制理論反推控制指令修正量。


本文采用方法1)預測終端狀態偏差,采用方法2)修正標稱控制量及標稱狀態量。其原因如下:對于終端偏差的預測而言,方法1)和方法2)的精度相當,且由于方法1)無需積分,計算速度更快。然而,方法1)成立的前提是基于一條準確的標稱軌跡,如果采用方法1)進行標稱軌跡的更新,將會導致標稱軌跡與實際軌跡的誤差不斷累計。因此,盡管方法2)涉及數值積分,速度較慢,但對于標稱軌跡的準確獲取仍是必要的。
(45)
(46)

(47)

將Φ展開,有
(48)

(49)

(50)
式中:wk為高斯型積分的權函數。

(51)
此二次規劃問題可利用成熟算法迭代求得數值解,此處給出其解析解的求解過程。定義拉格朗日函數L為
(52)

(53)
經代數運算,可得
(54)

(55)
P∈R6×6,設其第i行j列的元素為pij,則有
(56)


(57)


由于四元數±Q代表的是同一姿態變換,因此本文在進行線性偽譜模型預測控制修正時有以下特殊規定:規定姿態四元數Q的第1項q0≥0,若q0<0,取Q=-Q。作上述人為規定后,能夠避免MRPs在表示姿態時的奇異性,即當實際姿態偏差為小值時,MRPs參數偏差δσ也為小值,從而保證線性偽譜模型預測控制修正方法中的線性化前提成立。
線性偽譜模型預測能量最優姿態機動控制方法的具體實施步驟如下:
步驟1機動任務初始化:設置姿態機動的任務參數,包括姿態機動時間限制T,初始姿態角[γ0θ0ψ0],初始姿態角速度[ω10ω20ω30],終端姿態角[γfθfψf],終端姿態角速度[ω1fω2fω3f],控制幅值約束[U1maxU2maxU3max]等。


步驟4基準控制段:飛行器將按照標稱控制軌跡進行姿態機動控制,同時記錄當前的時間;若當前時間距上次更新檢查時間到達τcheck時,進入步驟5;若當前時間到達姿態機動限定時間T時,停止控制,完成姿態機動。


能量最優姿態機動控制仿真中使用的飛行器模型參數如表1所示。

表1 飛行器模型參數
姿態機動任務的初始姿態角為
[γ0θ0ψ0]=[-100 80 120](°)
姿態機動的初始角速度為
[ω10ω20ω30]=[0 0 0](°)/s
機動目標終端姿態角為
[γfθfψf]=[0 0 0](°)
機動目標終端角速度為
[ω1fω2fω3f]=[0 1 -0.5](°)/s
機動時間限制為
T=10 s
其中:γ為繞飛行器本體系OXbZ軸轉動的滾轉角;θ為繞飛行器本體系OYb軸轉動的俯仰角;ψ為繞飛行器本體系OZb軸轉動的偏航角。飛行器采用Z-Y-X的旋轉順序進行姿態變換。
取干擾力矩di~N(0,0.05),歐拉角測量誤差riangle~N(0,0.05)(實際中,從陀螺儀上獲取的姿態信息通常為歐拉角形式,故此處設定的測量誤差ri為歐拉角及角速度測量誤差,仿真中再通過參數變換得到MRPs的測量誤差),角速度測量誤差riangvel~N(0,0.01)。
取線性偽譜模型預測控制中的拉格朗日插值多項式的階數N=8。
指令更新間隔時間設為τcheck=2 s。
姿態角控制修正閾值設為0.01(單位:1),角速度控制修正閾值設為0.005(°)/s。
為對比控制效果,本文同時進行了4組的LQR姿態機動控制的仿真,分別為
1)Q=100×I6×6,R=0.01×I3×3,采用LQR方法跟蹤標稱軌跡。
2)Q=100×I6×6,R=0.01×I3×3,采用LQR方法直接進行機動。
3)Q=100×I6×6,R=0.1×I3×3,采用LQR方法跟蹤標稱軌跡。
4)Q=100×I6×6,R=0.1×I3×3,采用LQR方法直接進行機動。
本文采用的仿真環境為:Intel Core i5-4200M處理器,4 G內存,Windows 7 32位操作系統以及MATLAB R2013a。
設初始姿態角偏差為
[Δγ0Δθ0Δψ0]=[1.07 3.48 3.28](°)
初始角速度偏差為
[Δω10Δω20Δω30]=
[1.39 -1.83 1.75](°)/s
仿真結果如表2和圖2~圖4所示。
從仿真結果能夠看出,線性偽譜能量最優姿態機動控制能夠達到較高的控制精度,實現飛行器有限時間內的大角度姿態機動,并且相比于LQR控制消耗更少的能量,同時其控制指令更加平滑。
同時能夠看出,采用LQR跟蹤標稱軌跡的控制精度要優于直接控制。對于LQR控制而言,面臨著終端精度和消耗能量之間的權衡,若仿真參數R選取較大,盡管消耗能量變小,控制平滑,但終端精度下降;若仿真參數Q選取較大,終端狀態精度將提高,但消耗能量增大,同時控制量將發生振蕩。仿真表明,采用參數1的LQR跟蹤標稱軌跡的控制方法獲得的姿態機動終端精度與線性偽譜模型預測方法最為相近。
表2姿態機動單次仿真結果
Table2Singlesimulationresultsofattitudemaneuvers

控制方法[γf θf ψf]/(°)[ω1f ω2f ω3f]/((°)·s-1)Φ/(N2·m2·s)線性偽譜模型預測控制[0.15 -0.19 0.07][0.0081 0.94 -0.48]247.98LQR跟蹤標稱軌跡-參數1[0.26 -0.16 0.35][-0.097 1.04 -0.56]260.54LQR直接控制-參數1[7.76 17.7 -9.78][2.21 1.43 -3.28]489.22LQR跟蹤標稱軌跡-參數2[1.14 0.7 0.62][-0.13 0.75 -0.75]240.58LQR直接控制-參數2[8.7 23.4 -12.1][1.93 4.09 -2.76]249.75

圖2 能量最優姿態機動姿態角仿真曲線Fig.2 Simulation curves of attitude angles of fuel-optimal attitude maneuvers

圖3 能量最優姿態機動角速度仿真曲線Fig.3 Simulation curves of angular velocities of fuel-optimal attitude maneuvers

圖4 能量最優姿態機動控制力矩仿真曲線Fig.4 Simulation curves of control moment of fuel-optimal attitude maneuvers
通過蒙特卡羅仿真,可以統計與比較這幾種方法在計算精度、計算時間、消耗能量之間的差別。
為了在終端精度相當的前提下比較能量消耗,此處采用LQR跟蹤標稱軌跡(參數1)的方法與線性偽譜模型預測控制方法作蒙特卡羅仿真對比。設初始姿態角偏差滿足正態分布N(0,10),初始角速度偏差滿足正態分布N(0,5),做200次蒙特卡羅仿真,結果如表3和表4、圖5~圖8所示。

表3 姿態機動蒙特卡羅仿真終端精度

表4 姿態機動蒙特卡羅仿真仿真時間

圖5 姿態機動蒙特卡羅仿真姿態角曲線Fig.5 Curves of attitude angle maneuvers using Monte Carlo simulation

圖6 姿態機動蒙特卡羅仿真角速度曲線Fig.6 Curves of angular velocities of attitude maneuvers using Monte Carlo simulation

圖7 姿態機動蒙特卡羅仿真終端精度散布圖Fig.7 Scatter diagram of terminal accuracy of attitude maneuvers using Monte Carlo simulation

圖8 姿態機動蒙特卡羅仿真能量消耗對比圖Fig.8 Comparison of energy consumption of attitude maneuvers using Monte Carlo simulation
仿真中,線性偽譜能量最優姿態機動控制方法(節點數為8)進行一次指令修正更新的計算時間約為59 ms;LQR控制指令更新一次的計算時間約為14.6 ms。盡管線性偽譜控制指令修正的計算時間大于LQR控制,但該方法不需要每一時刻都計算更新指令,只有在預測偏差大于設定閾值的情況下才更新,在仿真中,指令更新數基本在1~3次的范圍內,故全過程仿真花費的時間兩者相差不大。因此,線性偽譜能量最優姿態機動控制方法能夠滿足飛行器上計算的實時性要求。
從仿真結果可以看出,線性偽譜能量最優姿態機動控制方法和LQR跟蹤方法都能滿足終端約束,能夠實現對初始偏差的修正,對干擾力矩以及測量誤差都有一定的抗干擾作用。然而,2種方法的姿態機動軌跡并不相同,LQR跟蹤方法趨向于盡快將當前狀態修正到標稱狀態,而線性偽譜能量最優姿態機動控制方法只關心終端狀態是否修正到目標狀態,對中間過程并不關心。實際上,當只關心終端約束時,并不需要立即將當前狀態修正到標稱狀態,而這種不必要的立即修正帶來了額外的能量消耗。
從控制的角度來看,線性偽譜能量最優姿態機動控制方法更具有優勢,因為其控制更加平穩,而LQR跟蹤的控制顯得更加振蕩。在考慮控制器的動態特性后,LQR的控制效果將會下降。
從圖8可以看出,線性偽譜能量最優姿態機動控制方法能量消耗小于LQR跟蹤控制,前者平均約能節省10%的能量。
本文基于線性偽譜模型預測控制,設計了線性偽譜能量最優姿態機動控制方法,該控制方法具有以下優勢:
1) 能夠實現限定時間內飛行器能量最優大角度姿態機動,相比于LQR跟蹤規劃軌跡的方法能夠節省約10%的能量消耗。
2) 該方法能夠得到的平滑且連續的控制量。
3) 該方法能夠快速地進行控制修正,滿足在線使用的要求。
4) 該方法的設計思路能推廣到更為一般的具有終端約束的微分動力學系統跟蹤問題上。
該方法目前還存在著以下缺陷:
1) 實際使用中,需根據飛行器的計算能力進行節點數、積分步長和控制精度之間的權衡。
2) 修正得到的控制量無法保證一定滿足過程約束,若直接在得到的修正控制量上增加限幅模塊,則得到的修正控制無法保證能量最優。同時可能導致額外的控制更新。
3) 本文提出的方法只適用于連續控制,不適用于Bang-Bang控制(開關式控制)。
第2)、3)點缺陷,將是未來的研究重點。