線性偽譜模型預測能量最優(yōu)姿態(tài)機動控制方法

馮逸駿， 陳萬春， 楊良

(北京航空航天大學 宇航學院, 北京 100083)

現(xiàn)代空間任務通常要求飛行器具有良好的姿態(tài)機動能力[1]?？臻g飛行器的姿態(tài)機動是指飛行器在外太空作大角度的姿態(tài)調(diào)整[2]。相比于人造衛(wèi)星、空間站等航天器，空間飛行器(如亞軌道飛行器、動能攔截器、軌道再入攻擊器等，以下簡稱飛行器)的姿態(tài)機動存在以下特點：①姿態(tài)機動角度大、時間短；②一般為主動控制，多采用反作用控制系統(tǒng)(Reaction Control Systems, RCS)，甚至是定推力的RCS進行控制；③關(guān)心能量消耗，在滿足精度條件下要求機動耗能盡可能??；④由于機動時間短，通常忽略如重力梯度、太陽光壓等擾動。以上特點為飛行器姿態(tài)機動控制系統(tǒng)的設計提出了新的挑戰(zhàn)。

飛行器在外太空姿態(tài)機動的過程中通?？煽醋鳛閯傮w。描述六自由度剛體姿態(tài)的系統(tǒng)通常是具有強耦合特性的多輸入多輸出(MIMO)非線性系統(tǒng)，因此飛行器能量最優(yōu)姿態(tài)機動控制問題實質(zhì)上是有限時域內(nèi)非線性系統(tǒng)最優(yōu)控制問題的一個特例。 針對此問題，Bharadwaj等提出了逆最優(yōu)控制(inverse optimal control)的方法[3]，通過構(gòu)造李雅普諾夫函數(shù)和求解HJB(Hamilton-Jacobi-Bellman)方程獲得反饋控制量，實現(xiàn)滿足性能指標最優(yōu)的動力學系統(tǒng)控制，并被廣泛用于飛行器能量最優(yōu)姿態(tài)機動和控制當中，如文獻[4-7]。然而，逆最優(yōu)控制方法滿足的是終端精度與能量消耗的加權(quán)和在無窮時域上的最優(yōu)，并非有限時間內(nèi)滿足終端條件的能量最優(yōu)控制。同時，逆最優(yōu)控制姿態(tài)機動的終端精度、機動時間和消耗能量取決于若干控制參數(shù)的調(diào)節(jié)，在實際使用中有所不便。

為此，本文提出了一種基于模型預測控制(Model Predictive Control,MPC)和線性偽譜的能量最優(yōu)姿態(tài)機動控制方法。

模型預測控制是一種基于滾動優(yōu)化的在線控制策略，具有對模型要求低、抗干擾性好、魯棒性強等優(yōu)點，能夠在優(yōu)化性能指標的同時處理各種約束條件，得到了工程技術(shù)人員和理論研究者的重視，并被應用到飛行器姿態(tài)控制問題的研究當中[1,8-9]。近年來，快速MPC逐漸引起關(guān)注。模型預測控制通常涉及到當前狀態(tài)線性化后的局部最優(yōu)控制問題，一般地，求解該最優(yōu)控制問題可轉(zhuǎn)化為求解一個兩點邊值問題(Two-Point Boundary Value Problem,TPBVP)，專家學者為提高控制指令求解速率和控制精度做了大量工作[10]，針對具有強終端約束以及二次性能指標的控制問題，提出了模型預測靜態(tài)規(guī)劃(Model Predictive Static Programming,MPSP)非線性最優(yōu)控制方法[11]，并將其成功運用在飛行器再入制導領域[12]。盡管MPSP方法能獲得全局最優(yōu)解，但其在計算過程中需要采用大量離散節(jié)點以保證積分精度，導致計算效率較低。針對該問題，Yang等[13]提出了線性高斯偽譜模型預測控制(Linear Gauss Pseudo-spectral Model Predictive Control,LGPMPC)方法，該方法綜合了非線性近似模型預測控制、線性二次最優(yōu)控制以及高斯偽譜法，采用較MPSP更少的離散節(jié)點達到更高的計算效率及精度，被應用于大氣層外制導[14]。針對多段問題，Yang等[15]又在LGPMPC的基礎上提出了多段線性偽譜模型預測控制(Multi-segment Linear Pseudo-spectral Model Predictive Control,MLPMPC)方法，并被成功應用于飛行器再入段制導。

本文在LGPMPC和MLPMPC方法的基礎上，針對飛行器姿態(tài)機動問題，利用KKT (Karush-Kuhn-Tucher)條件和線性偽譜推導了能量最優(yōu)控制指令修正量的解析表達式，從而得到基于線性偽譜模型預測控制的能量最優(yōu)姿態(tài)機動控制方法。該方法能夠同時滿足多個終端約束，并保證修正后的控制指令仍然滿足能量最優(yōu)。基于線性偽譜模型預測控制的能量最優(yōu)姿態(tài)機動控制方法的基本策略為：根據(jù)飛行任務離線優(yōu)化出初始姿態(tài)機動軌跡，實際飛行過程中在線修正標稱控制以保證姿態(tài)機動終端精度，同時迭代更新姿態(tài)機動標稱軌跡。數(shù)值仿真結(jié)果表明，該方法能夠在限定時間內(nèi)實現(xiàn)能量最優(yōu)大角度姿態(tài)機動，在近似精度的前提下，相比于傳統(tǒng)線性二次型調(diào)節(jié)器(Linear Quadratic Regulator,LQR)跟蹤節(jié)省約10%的能量消耗。

1 飛行器模型及姿態(tài)機動問題

本文采用修正羅德里格斯參數(shù)(Modified Rodrigues Parameters,MRPs)描述飛行器姿態(tài)?？紤]一般的剛體飛行器，其姿態(tài)運動學和動力學方程可描述為

(1)

(2)

G(σ)∈R3×3為飛行器姿態(tài)運動學矩陣，定義為

(3)

式中：I3×3為3×3單位矩陣。飛行器姿態(tài)機動過程中，控制力矩存在如式(4)幅值約束：

|ui|≤Uimaxi=1,2,3

(4)

式中：Uimax>0為每個姿態(tài)控制通道上的最大力矩幅值。

設狀態(tài)量x=[σ1σ2σ3ω1ω2ω3]T，則飛行器姿態(tài)運動系統(tǒng)可采用狀態(tài)空間的形式描述為

(5)

(6)

式中：t0和tf分別為飛行器姿態(tài)機動的起始時刻和終端時刻；φ為末值型項指標；φ為積分型項指標。此處性能指標Φ只考慮能量最優(yōu)，故只與控制量相關(guān)，一般地，取

(7)

2 線性偽譜能量最優(yōu)姿態(tài)機動控制

線性偽譜模型預測控制方法的主要思想為：將非線性動力學方程在標稱狀態(tài)量附近進行擬線性化，建立一個以偏差為自變量的線性微分方程，通過拉格朗日插值多項式對該線性微分方程的狀態(tài)量進行逼近，將微分動力學約束通過正交配點轉(zhuǎn)化為一組多變量代數(shù)約束，最終將非線性方程的積分問題轉(zhuǎn)化為一個連續(xù)求解線性代數(shù)方程組的問題，從而得到非線性動力學過程的終端狀態(tài)與狀態(tài)偏差、控制偏差的等式關(guān)系，以實現(xiàn)偏差修正[13-14]。

基于線性偽譜模型預測控制的能量最優(yōu)姿態(tài)機動控制策略需要基于一條標稱機動軌跡，因此須先求解最優(yōu)控制問題，規(guī)劃得到標稱機動軌跡。本文采用高斯偽譜法將能量最優(yōu)姿態(tài)機動問題離散為非線性規(guī)劃問題，采用SNOPT工具包求解該非線性規(guī)劃問題。標稱機動軌跡的具體求解及優(yōu)化過程，并非本文研究重點，此處不再贅述。

線性偽譜模型預測能量最優(yōu)姿態(tài)機動控制方法可以分為終端狀態(tài)解析預測和控制指令修正兩部分。

2.1 線性偽譜姿態(tài)機動終端狀態(tài)解析預測

一般地，考慮具有終端約束的姿態(tài)機動非線性動力學方程如下：

(8)

(9)

(10)

(11)

式中：A(t)∈R6×6為狀態(tài)誤差傳播矩陣，B(t)∈R6×3為控制誤差傳播矩陣。對式(1)～式(3)描述的飛行器姿態(tài)動力學系統(tǒng)而言，狀態(tài)誤差傳播矩陣和控制誤差傳播矩陣的具體形式為

(12)

其中：A11∈R3×3，A12∈R3×3，A22∈R3×3，具體表達式分別為

(13)

(14)

(15)

式中：K、S1、S2和S3為中間變量，其表達式分別為

(16)

(17)

(18)

(19)

對狀態(tài)量和控制量進行偽譜離散時，通常有3種高斯正交節(jié)點可供選擇：LG(Legendre Gauss)節(jié)點、LGR(Legendre Gauss Radau)節(jié)點、LGL(Legendre Gauss Lobatto)節(jié)點。其中，LG節(jié)點不含終端點，LGR節(jié)點含有一個終端點(通常為左端歸一化時間τ=-1)，LGL節(jié)點含有2個終端點，如圖1所示?？紤]到工程實際中的控制連續(xù)性，此處選取LGR節(jié)點對狀態(tài)量和控制量進行離散。

首先，將實際機動時間t∈[t0,tf]映射到歸一化時間τ∈[-1,1]上：

(20)

則歸一化后的誤差傳播動力學方程為

(21)

式中：p為歸一化轉(zhuǎn)換變量，其表達式為

(22)

(23)

同樣地，本文得到控制量的插值擬合：

(24)

拉格朗日插值多項式滿足以下性質(zhì)：

(25)

故有

(26)

通過對狀態(tài)量求導可得

(27)

式中：微分逼近矩陣D∈RN×(N+1)是通過對拉格朗日插值多項式的各個元素分別求導獲得的，其具體表達式為

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

(28)

(29)

式中：k=1,2,…,N。

(30)

重新組合微分逼近矩陣D，可得到預測狀態(tài)偏差序列和控制修正序列表示的關(guān)系式如下：

(31)

其中:

(32)

(33)

(34)

矩陣A*和B*的表達式分別為

(35)

(36)

則其他LGR節(jié)點上的各個狀態(tài)量可表示為

(37)

(38)

(39)

其中：Wx為高斯型積分公式的權(quán)函數(shù)矩陣，可表示為

(40)

(41)

(42)

式中：Kx∈R6×6，Ku∈R6×3N，具體表達式分別為

Kx=w0pA(τ0)-

pA′]-1D1

(43)

pA′]-1B*

(44)

至此，本文獲得了預測終端狀態(tài)偏差關(guān)于初始狀態(tài)偏差和控制修正的解析等式關(guān)系。

2.2 線性偽譜能量最優(yōu)姿態(tài)機動控制指令修正

在求得預測終端狀態(tài)偏差關(guān)于初始狀態(tài)偏差和控制修正的解析等式關(guān)系后，可以根據(jù)最優(yōu)控制理論反推控制指令修正量。

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

(45)

(46)

(47)

將Φ展開，有

(48)

(49)

(50)

式中：wk為高斯型積分的權(quán)函數(shù)。

(51)

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

(52)

(53)

經(jīng)代數(shù)運算，可得

(54)

(55)

P∈R6×6，設其第i行j列的元素為pij，則有

(56)

(57)

由于四元數(shù)±Q代表的是同一姿態(tài)變換，因此本文在進行線性偽譜模型預測控制修正時有以下特殊規(guī)定：規(guī)定姿態(tài)四元數(shù)Q的第1項q0≥0，若q0<0，取Q=-Q。作上述人為規(guī)定后，能夠避免MRPs在表示姿態(tài)時的奇異性，即當實際姿態(tài)偏差為小值時，MRPs參數(shù)偏差δσ也為小值，從而保證線性偽譜模型預測控制修正方法中的線性化前提成立。

2.3 能量最優(yōu)姿態(tài)機動控制方法實施步驟

線性偽譜模型預測能量最優(yōu)姿態(tài)機動控制方法的具體實施步驟如下：

步驟1機動任務初始化：設置姿態(tài)機動的任務參數(shù)，包括姿態(tài)機動時間限制T，初始姿態(tài)角[γ0θ0ψ0]，初始姿態(tài)角速度[ω10ω20ω30]，終端姿態(tài)角[γfθfψf]，終端姿態(tài)角速度[ω1fω2fω3f]，控制幅值約束[U1maxU2maxU3max]等。

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

3 能量最優(yōu)姿態(tài)機動控制仿真

3.1 飛行器模型參數(shù)與控制參數(shù)設定

能量最優(yōu)姿態(tài)機動控制仿真中使用的飛行器模型參數(shù)如表1所示。

表1 飛行器模型參數(shù)

姿態(tài)機動任務的初始姿態(tài)角為

[γ0θ0ψ0]=[-100 80 120](°)

姿態(tài)機動的初始角速度為

[ω10ω20ω30]=[0 0 0](°)/s

機動目標終端姿態(tài)角為

[γfθfψf]=[0 0 0](°)

機動目標終端角速度為

[ω1fω2fω3f]=[0 1 -0.5](°)/s

機動時間限制為

T=10 s

其中：γ為繞飛行器本體系OXbZ軸轉(zhuǎn)動的滾轉(zhuǎn)角；θ為繞飛行器本體系OYb軸轉(zhuǎn)動的俯仰角；ψ為繞飛行器本體系OZb軸轉(zhuǎn)動的偏航角。飛行器采用Z-Y-X的旋轉(zhuǎn)順序進行姿態(tài)變換。

取干擾力矩di～N(0,0.05)，歐拉角測量誤差riangle～N(0,0.05)(實際中，從陀螺儀上獲取的姿態(tài)信息通常為歐拉角形式，故此處設定的測量誤差ri為歐拉角及角速度測量誤差，仿真中再通過參數(shù)變換得到MRPs的測量誤差)，角速度測量誤差riangvel～N(0,0.01)。

取線性偽譜模型預測控制中的拉格朗日插值多項式的階數(shù)N=8。

指令更新間隔時間設為τcheck=2 s。

姿態(tài)角控制修正閾值設為0.01(單位：1)，角速度控制修正閾值設為0.005(°)/s。

為對比控制效果，本文同時進行了4組的LQR姿態(tài)機動控制的仿真，分別為

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方法直接進行機動。

本文采用的仿真環(huán)境為：Intel Core i5-4200M處理器，4 G內(nèi)存，Windows 7 32位操作系統(tǒng)以及MATLAB R2013a。

3.2 能量最優(yōu)姿態(tài)機動單次仿真

設初始姿態(tài)角偏差為

[Δγ0Δθ0Δψ0]=[1.07 3.48 3.28](°)

初始角速度偏差為

[Δω10Δω20Δω30]=

[1.39 -1.83 1.75](°)/s

仿真結(jié)果如表2和圖2～圖4所示。

從仿真結(jié)果能夠看出，線性偽譜能量最優(yōu)姿態(tài)機動控制能夠達到較高的控制精度，實現(xiàn)飛行器有限時間內(nèi)的大角度姿態(tài)機動，并且相比于LQR控制消耗更少的能量，同時其控制指令更加平滑。

同時能夠看出，采用LQR跟蹤標稱軌跡的控制精度要優(yōu)于直接控制。對于LQR控制而言，面臨著終端精度和消耗能量之間的權(quán)衡，若仿真參數(shù)R選取較大，盡管消耗能量變小，控制平滑，但終端精度下降；若仿真參數(shù)Q選取較大，終端狀態(tài)精度將提高，但消耗能量增大，同時控制量將發(fā)生振蕩。仿真表明，采用參數(shù)1的LQR跟蹤標稱軌跡的控制方法獲得的姿態(tài)機動終端精度與線性偽譜模型預測方法最為相近。

表2姿態(tài)機動單次仿真結(jié)果

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跟蹤標稱軌跡-參數(shù)1[0.26 -0.16 0.35][-0.097 1.04 -0.56]260.54LQR直接控制-參數(shù)1[7.76 17.7 -9.78][2.21 1.43 -3.28]489.22LQR跟蹤標稱軌跡-參數(shù)2[1.14 0.7 0.62][-0.13 0.75 -0.75]240.58LQR直接控制-參數(shù)2[8.7 23.4 -12.1][1.93 4.09 -2.76]249.75

圖2 能量最優(yōu)姿態(tài)機動姿態(tài)角仿真曲線Fig.2 Simulation curves of attitude angles of fuel-optimal attitude maneuvers

圖3 能量最優(yōu)姿態(tài)機動角速度仿真曲線Fig.3 Simulation curves of angular velocities of fuel-optimal attitude maneuvers

圖4 能量最優(yōu)姿態(tài)機動控制力矩仿真曲線Fig.4 Simulation curves of control moment of fuel-optimal attitude maneuvers

通過蒙特卡羅仿真，可以統(tǒng)計與比較這幾種方法在計算精度、計算時間、消耗能量之間的差別。

3.3 能量最優(yōu)姿態(tài)機動蒙特卡羅仿真

為了在終端精度相當?shù)那疤嵯卤容^能量消耗，此處采用LQR跟蹤標稱軌跡(參數(shù)1)的方法與線性偽譜模型預測控制方法作蒙特卡羅仿真對比。設初始姿態(tài)角偏差滿足正態(tài)分布N(0,10)，初始角速度偏差滿足正態(tài)分布N(0,5)，做200次蒙特卡羅仿真，結(jié)果如表3和表4、圖5～圖8所示。

表3 姿態(tài)機動蒙特卡羅仿真終端精度

表4 姿態(tài)機動蒙特卡羅仿真仿真時間

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

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

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

圖8 姿態(tài)機動蒙特卡羅仿真能量消耗對比圖Fig.8 Comparison of energy consumption of attitude maneuvers using Monte Carlo simulation

仿真中，線性偽譜能量最優(yōu)姿態(tài)機動控制方法(節(jié)點數(shù)為8)進行一次指令修正更新的計算時間約為59 ms；LQR控制指令更新一次的計算時間約為14.6 ms。盡管線性偽譜控制指令修正的計算時間大于LQR控制，但該方法不需要每一時刻都計算更新指令，只有在預測偏差大于設定閾值的情況下才更新，在仿真中，指令更新數(shù)基本在1～3次的范圍內(nèi)，故全過程仿真花費的時間兩者相差不大。因此，線性偽譜能量最優(yōu)姿態(tài)機動控制方法能夠滿足飛行器上計算的實時性要求。

從仿真結(jié)果可以看出，線性偽譜能量最優(yōu)姿態(tài)機動控制方法和LQR跟蹤方法都能滿足終端約束，能夠?qū)崿F(xiàn)對初始偏差的修正，對干擾力矩以及測量誤差都有一定的抗干擾作用。然而，2種方法的姿態(tài)機動軌跡并不相同，LQR跟蹤方法趨向于盡快將當前狀態(tài)修正到標稱狀態(tài)，而線性偽譜能量最優(yōu)姿態(tài)機動控制方法只關(guān)心終端狀態(tài)是否修正到目標狀態(tài)，對中間過程并不關(guān)心。實際上，當只關(guān)心終端約束時，并不需要立即將當前狀態(tài)修正到標稱狀態(tài)，而這種不必要的立即修正帶來了額外的能量消耗。

從控制的角度來看，線性偽譜能量最優(yōu)姿態(tài)機動控制方法更具有優(yōu)勢，因為其控制更加平穩(wěn)，而LQR跟蹤的控制顯得更加振蕩。在考慮控制器的動態(tài)特性后，LQR的控制效果將會下降。

從圖8可以看出，線性偽譜能量最優(yōu)姿態(tài)機動控制方法能量消耗小于LQR跟蹤控制，前者平均約能節(jié)省10%的能量。

4 結(jié) 論

本文基于線性偽譜模型預測控制，設計了線性偽譜能量最優(yōu)姿態(tài)機動控制方法，該控制方法具有以下優(yōu)勢：

1) 能夠?qū)崿F(xiàn)限定時間內(nèi)飛行器能量最優(yōu)大角度姿態(tài)機動，相比于LQR跟蹤規(guī)劃軌跡的方法能夠節(jié)省約10%的能量消耗。

2) 該方法能夠得到的平滑且連續(xù)的控制量。

3) 該方法能夠快速地進行控制修正，滿足在線使用的要求。

4) 該方法的設計思路能推廣到更為一般的具有終端約束的微分動力學系統(tǒng)跟蹤問題上。

該方法目前還存在著以下缺陷：

1) 實際使用中，需根據(jù)飛行器的計算能力進行節(jié)點數(shù)、積分步長和控制精度之間的權(quán)衡。

2) 修正得到的控制量無法保證一定滿足過程約束，若直接在得到的修正控制量上增加限幅模塊，則得到的修正控制無法保證能量最優(yōu)。同時可能導致額外的控制更新。

3) 本文提出的方法只適用于連續(xù)控制，不適用于Bang-Bang控制(開關(guān)式控制)。

第2)、3)點缺陷，將是未來的研究重點。