石青鑫 戈新生

(北京信息科技大學(xué)機電工程學(xué)院, 北京 100192)

引言

欠驅(qū)動航天器是指僅依靠兩個執(zhí)行機構(gòu)控制實現(xiàn)姿態(tài)機動的航天器.Crouch[1]研究了欠驅(qū)動航天器的可控性,證明了在系統(tǒng)的動量矩等于零時,欠驅(qū)動航天器可控.Morin[2]設(shè)計了指數(shù)收斂的時變狀態(tài)反饋穩(wěn)定控制器,實現(xiàn)了欠驅(qū)動航天器的姿態(tài)穩(wěn)定控制.戈新生等[3]利用Ritz理論和粒子群算法,研究了帶有兩動量輪的欠驅(qū)動航天器和太陽帆板展開過程的姿態(tài)控制問題.王冬霞等[4]利用分層滑模控制方法,為帶兩個控制執(zhí)行機構(gòu)的欠驅(qū)動剛體航天器的姿態(tài)控制系統(tǒng)設(shè)計了一種三軸穩(wěn)定控制器.莊宇飛等[5]應(yīng)用偽譜法解決了欠驅(qū)動剛性航天器的時間最優(yōu)軌跡規(guī)劃問題,表明偽譜法能夠較好地滿足各種約束條件,而且計算精度高、速度快,具有良好的實時性.近年來,偽譜法越來越成為解決最優(yōu)控制問題(OCP)和非線性方程的重要方法,其根本原因是NLP問題的Karush-Kuhn-Tucker(KKT)條件與離散的哈密頓邊值問題的一階最優(yōu)條件等價得到了證明[6].偽譜法主要包括Legendre偽譜法、Gauss偽譜法、Radau偽譜法、Chebyshev偽譜法,其中Legendre偽譜法和Gauss偽譜法使用較為普遍.唐小軍等[7]提出以Chebyshev-Gauss(CG)偽譜法解決OCP的直接軌跡優(yōu)化,嚴(yán)格推導(dǎo)了約束乘子估計與NLP問題的KKT條件的等價性.

航天器的姿態(tài)機動可以通過噴氣推力器或動量輪實現(xiàn).使用CG偽譜法研究欠驅(qū)動航天器姿態(tài)最優(yōu)控制問題,采用CG配置點,利用 Clenshaw-Curtis積分近似得到性能指標(biāo)函數(shù)中的積分項,并用快速傅里葉變換(FFT)以快速有效計算相應(yīng)的積分權(quán)；為提高計算效率和數(shù)值穩(wěn)定性,應(yīng)用重心拉格朗日插值代替經(jīng)典的拉格朗日插值去逼近狀態(tài)變量和控制變量,并采用三角函數(shù)來減小計算配點數(shù)較大時的狀態(tài)微分矩陣帶來的誤差,將連續(xù)OCP離散為具有代數(shù)約束的參數(shù)優(yōu)化NLP問題,再通過SQP算法求解得到控制輸入的規(guī)律.通過數(shù)值仿真得出,對兩類欠驅(qū)動航天器的姿態(tài)最優(yōu)控制均能達(dá)到設(shè)計的零邊界控制要求,得到的姿態(tài)最優(yōu)曲線與反向驗證得到的曲線幾乎完全重疊.

1 欠驅(qū)動航天器模型

設(shè)航天器系統(tǒng)由一個航天器主體B及兩個動量輪Wi(i=1,2)組成,并且兩輪的自轉(zhuǎn)軸分別與航天器主體的前兩個主軸重合(如圖1所示),動量輪質(zhì)心為Ci(i=1,2).以航天器系統(tǒng)質(zhì)心建立慣性系(O-XYZ),以動量輪的質(zhì)心建立連體坐標(biāo)系分別為(Oi-xiyizi)(i=1,2),并設(shè)x1,y2分別為兩動量輪的自轉(zhuǎn)軸.

圖1 帶兩動量輪的欠驅(qū)動航天器簡化圖Fig.1 Simplified representation of underactuated spacecraft with two momentum wheels

使用歐拉角ψ,θ,φ分別代表偏航向角,俯仰角,滾動角,以3-2-1次序旋轉(zhuǎn),得到方向矩陣[8]

(1)

式中sin(*)=s*,cos(*)=c*,航天器姿態(tài)運動學(xué)方程為[9,10]：

(2)

(3)

將式(2)展開可得以歐拉角描述姿態(tài)的運動學(xué)方程：

(4)

設(shè)航天器系統(tǒng)對質(zhì)心O的總動量矩為H,可表示為：

(5)

其中JB=diag(JB1,JB2,JB3)和JWi分別為航天器主體、兩動量輪對質(zhì)心O轉(zhuǎn)動慣量矩陣,JCi為動量輪對其質(zhì)心的轉(zhuǎn)動慣量矩陣,Ωi為動量輪相對航天器主體轉(zhuǎn)動角速度矢量.

由于動量輪相對航天器主體作定軸轉(zhuǎn)動,式(5)右側(cè)第三項可簡化為：

(6)

(7)

其中J為航天器系統(tǒng)對質(zhì)心O的轉(zhuǎn)動慣量矩陣,表示如下：

(8)

根據(jù)動量矩定理有：

(9)

式中M表示施加在系統(tǒng)質(zhì)心的外力矩,將式(7)對時間求導(dǎo)：

(10)

(11)

當(dāng)航天器系統(tǒng)以動量輪作為執(zhí)行機構(gòu)時,式(11)中M=0,此時令動量輪的角加速度作為系統(tǒng)的控制變量,即：

(12)

則帶動量輪的欠驅(qū)動航天器系統(tǒng)動力學(xué)模型為：

(13)

其中：

(14)

該系統(tǒng)對應(yīng)的狀態(tài)變量為:

(15)

該系統(tǒng)對應(yīng)的狀態(tài)變量為:

2 欠驅(qū)動航天器姿態(tài)最優(yōu)控制問題

為使欠驅(qū)動航天器既實現(xiàn)期望的姿態(tài)位置,又滿足執(zhí)行工作任務(wù)的其它要求(譬如運行時間最短、消耗燃料最少等),最優(yōu)控制問題由此而生,在此選取Bolza形式的最優(yōu)控制：

(16)

上述方程依次代表：性能指標(biāo)函數(shù)、狀態(tài)約束方程、邊界約束方程和路徑約束方程.其中x∈n為狀態(tài)變量,用u∈m代表上文提到的控制變量和M,t為實際任意時間,t0為實際初始時間,tf為實際終端時間.函數(shù)Φ和g是標(biāo)量,f∈n為n維向量函數(shù),φ∈q為q維向量函數(shù),C∈c為c維向量函數(shù).

根據(jù)可控性秩條件(Chow定理)[11]可知：控制系統(tǒng)存在一個滿秩可控性李代數(shù),原則上能求解系統(tǒng)的姿態(tài)機動問題.設(shè)控制系統(tǒng)存在最優(yōu)解u*∈L2([0,T]),其中L2([0,T])為可測向量函數(shù)u(t),t∈[0,T]構(gòu)成的Hilbert空間.根據(jù)最優(yōu)控制原理,選取優(yōu)化指標(biāo)為：

(17)

對于帶動量輪的航天器,給定系統(tǒng)初末位形、角速度和動量輪的角速度,即x0,xf∈8；而對于帶推力器的航天器,給定初末位形和角速度,即x0,xf∈6.通過目標(biāo)函數(shù)尋求控制輸入u(t)∈2,t∈[0,T]T,從而可確定系統(tǒng)從x0機動到xf的姿態(tài)運動軌跡.控制航天器的執(zhí)行機構(gòu)在啟動和關(guān)閉瞬間的控制值應(yīng)均為零,即u0=uf=0.另外,控制值不可能無限大,應(yīng)設(shè)置最大值,因此還需考慮

3 Chebyshev-Gauss偽譜法

CG偽譜法通過一系列變換將連續(xù)最優(yōu)控制問題離散為具有代數(shù)約束的NLP問題.

首先進(jìn)行時域變換,需要引入新的時間變量τ∈[-1,1],將時間區(qū)間[t0,tf]轉(zhuǎn)換為[-1,1],定義時間變量t為：

(18)

其次對狀態(tài)變量與控制變量作如下處理.取K個CG點以及τ0=-1作為CG偽譜法的節(jié)點,其中CG點[12]的定義為：

(19)

CG點在區(qū)間(-1,1)上的分布特點為兩端稠密,中間稀疏,能有效抑制拉格朗日插值時兩端容易出現(xiàn)的Runge現(xiàn)象.用此K+1個節(jié)點構(gòu)造拉格朗日插值多項式,并以此為基函數(shù)來逼近狀態(tài)變量有：

(20)

為提高計算效率和數(shù)值穩(wěn)定性,采用重心拉格朗日插值[13]：

(21)

其中重心權(quán)ξi的定義為：

(22)

對于式(22),在數(shù)值計算τi-τj時會產(chǎn)生浮點誤差[13],為避免誤差作如下處理,設(shè)：

(23)

考慮到τK+1=1,由式(22)和(23)可得：

ξi=(τi-τK+1)ξi′

(24)

(25)

基于時域變換、狀態(tài)變量與控制變量的近似變換,式(16)寫作：

(26)

然后將運動學(xué)與動力學(xué)微分方程進(jìn)行離散,在CG偽譜法中,由全局正交多項式近似狀態(tài)變量,可通過對式(20)求導(dǎo)來近似得到：

(27)

式中τk為CG點.文中采用重心拉格朗日插值,因此微分矩陣D∈RK×(K+1)表示如下,該矩陣具有較好的數(shù)值穩(wěn)定性[14]：

(28)

式中k=1,…,K,i=0,…,K.聯(lián)立式(27)、(28)以及(26)中的狀態(tài)約束方程,得到離散方程：

(29)

再將終端狀態(tài)約束離散化,由于式(26)中包含了終端狀態(tài)約束,而式(20)的時間區(qū)間為[-1,1),因而狀態(tài)變量的離散表達(dá)式中未包含末端時間節(jié)點X(τf),對(16)中的狀態(tài)約束方程在區(qū)間[-1,1]上積分：

(30)

將式(30)左邊展開得：

(31)

對式(31)中的積分項采用Clenshaw-Curtis積分[15]求其近似值：

(32)

(33)

(34)

(35)

同樣利用Clenshaw-Curtis積分方法將指標(biāo)函數(shù)離散化,近似式(26)中的性能指標(biāo)函數(shù),可得：

J= Φ(X0,t0,Xf,tf)+

(36)

考慮到系統(tǒng)模型有兩個控制變量,即U1K,U2K∈K,再根據(jù)式(17),式(36)變?yōu)椋?/p>

(37)

其中μcc為Clenshaw-Curtis積分權(quán)向量.設(shè)置狀態(tài)變量和控制變量初末值及控制約束條件：

x0-X0=0,xf-Xf=0;

(38)

式中X0,Xf分別代表近似狀態(tài)變量的初末值,相應(yīng)的控制變量初末值分別為U0,Uf.Umax為控制變量的最大值.

至此,基于CG偽譜法的最優(yōu)控制問題轉(zhuǎn)化為NLP問題,即在滿足式(29)、(32)和式(38)的約束條件下尋求最優(yōu)控制輸入以使式(37)取得最小值.

4 數(shù)值仿真

取系統(tǒng)質(zhì)量幾何參數(shù)為[17]：

JB=diag(86.215,85.07,113.565)kg·m2

JW1=diag(0.5,0.45,0.45)kg·m2

JW2=diag(0.45,0.5,0.45)kg·m2

J1=J2=0.5kg·m2

假設(shè)欠驅(qū)動航天器系統(tǒng)進(jìn)行rest-to-rest的姿態(tài)機動,設(shè)置機動時間為20s,給定航天器姿態(tài)角ψ,θ,φ的初始值依次為0,0,0,末端值為0,π/6,0,控制輸入初始和終端值均為零.優(yōu)化算法采用從可行解到最優(yōu)解的串行優(yōu)化策略,可行解選取CG點個數(shù)為6個,最優(yōu)解選取CG點個數(shù)為60個,根據(jù)上文的理論公式,在內(nèi)存為4.0G,CPU頻率為2.1GHz的Windows7操作系統(tǒng)中使用MATLAB R2015b編程仿真.

帶動量輪的欠驅(qū)動航天器姿態(tài)最優(yōu)控制結(jié)果如圖2～圖5所示.

圖2 帶動量輪欠驅(qū)動航天器姿態(tài)角的機動曲線Fig.2 Optimal trajectory of the underactuated spacecraft with momentum wheels

圖中CGPM代表CG偽譜法.為驗證文中姿態(tài)優(yōu)化曲線的合理性及有效性,將求解得到的圖5中的控制變量代入到狀態(tài)方程式(4)和式(13)中使用ode45進(jìn)行積分,積分后得出的狀態(tài)變量結(jié)果如圖2～圖4中的藍(lán)色虛線所示,可以看出得到的積分曲線與用CG偽譜法優(yōu)化得到的曲線幾乎是重疊的,以圖2姿態(tài)角ψ曲線圖為例,觀察其局部放大圖可以發(fā)現(xiàn)兩曲線誤差在10-2數(shù)量級,可以看出該算法求解得到的姿態(tài)運動以及控制輸入都是合理且有效的.

圖3 帶動量輪欠驅(qū)動航天器姿態(tài)角速度的機動曲線Fig.3 Optimal angular velocity of the underactuated spacecraft with momentum wheels

為直觀體現(xiàn)CG偽譜法的優(yōu)勢,文中同時給出采用Gauss偽譜法計算的結(jié)果.Gauss偽譜法的配置點為Legendre-Gauss(LG)點；采用Gauss積分來離散性能指標(biāo)函數(shù)中的積分項；應(yīng)用經(jīng)典拉格朗日插值去逼近狀態(tài)和控制變量.Gauss偽譜法同樣采用串行優(yōu)化策略,可行解選取LG點個數(shù)為6個,最優(yōu)解選取LG點個數(shù)為60個,狀態(tài)和控制變量的參數(shù)設(shè)置均與CG偽譜法的相同.CG偽譜法和Gauss偽譜法可行解的初始猜測值均采用MATLAB的隨機函數(shù)生成,結(jié)果如表1所示.表中的參數(shù)Itf和Jf分別表示可行解求解階段的迭代次數(shù)以及目標(biāo)函數(shù)值；參數(shù)Ito和Jo則分別代表最優(yōu)解求解階段的迭代次數(shù)以及目標(biāo)函數(shù)值；而參數(shù)Ae和tCPU則分別表示最終的末端姿態(tài)最大誤差以及CPU總的運行時間.表1直觀地表明,CG偽譜法的目標(biāo)函數(shù)的可行解和最優(yōu)解值均比Gauss偽譜法的小,求解時間比Gauss偽譜法明顯縮短.

圖4 兩動量輪轉(zhuǎn)動角速度曲線Fig.4 Optimal angular velocity of two momentum wheels

圖5 最優(yōu)控制輸入(兩動量輪角加速度)曲線Fig.5 Optimal control input for two momentum wheel accelerations

表1 CG偽譜法與Gauss偽譜法的比較Table 1 Comparisonof CG Pseudospectral Method and Gauss Pseudospectral Method

帶推力器的欠驅(qū)動航天器姿態(tài)最優(yōu)控制結(jié)果如圖6～圖8所示.

圖6 帶推力器欠驅(qū)動航天器姿態(tài)角的機動曲線Fig.6 Optimal trajectory of the underactuated spacecraft with thrusters

圖7 帶推力器欠驅(qū)動航天器姿態(tài)角速度的機動曲線Fig.7 Optimal angular velocity of the underactuated spacecraft with thrusters

從圖6～圖8中的曲線可看到,姿態(tài)角和姿態(tài)角速度均可精確機動到預(yù)定的末端值,且相應(yīng)的軌跡曲線是光滑的.

5 結(jié)論

本文利用CG偽譜法解決了帶動量輪、推力器兩類欠驅(qū)動航天器姿態(tài)機動最優(yōu)控制問題.基于CG偽譜法把系統(tǒng)的連續(xù)姿態(tài)運動規(guī)劃轉(zhuǎn)化成NLP問題,再用求解非線性約束優(yōu)化問題的SQP算法運算.通過對這兩類模型的數(shù)值仿真,得到了系統(tǒng)從初始姿態(tài)機動到終端姿態(tài)的運動軌跡,終端角速度也達(dá)到了預(yù)定值,而且最優(yōu)控制輸入的初末值均等于零,可以說明該方法的可行性.將求得的最優(yōu)控制解代回系統(tǒng)的姿態(tài)運動方程中,利用數(shù)值積分求解狀態(tài)變量的結(jié)果,可以發(fā)現(xiàn)兩種方法求解得到的曲線幾乎完全重疊,進(jìn)而反向驗證了CG偽譜法對兩類系統(tǒng)的姿態(tài)最優(yōu)控制的合理性和有效性.另外,使用 Gauss偽譜法進(jìn)行仿真對比,看出CG偽譜法的計算效率明顯更高.

圖8 推力器的最優(yōu)控制輸入(力矩)曲線Fig.8 Optimal control input for the thrusters

1Crouch P E. Application of geometric control theory to rigid body models.IEEETransactionsonAutomaticControl, 1984,29(4):87～95

2Morin P,Samson C. Time-varying exponential stabilization of a rigid spacecraft with two control torques.IEEETransactionsonAutomaticControl, 1997,42(4):528～534

3戈新生,陳立群,劉延柱. 帶有兩個動量飛輪剛體航天器的姿態(tài)非完整運動規(guī)劃問題. 控制理論與應(yīng)用, 2004,21(5):781～784 (Ge X S, Chen L Q, Liu Y Z. Nonholonomic motion planning for the attitude of rigid spacecraft with two momentum wheel actuator.ControlTheory&Applications, 2004,21(5):781～784 (in Chinese))

4王冬霞,賈英宏,金磊. 欠驅(qū)動航天器姿態(tài)穩(wěn)定的分層滑?？刂破髟O(shè)計. 宇航學(xué)報, 2013,34(1):17～24 (Wang D X, Jia Y H, Jin L. Hierarchical sliding-mode control for attitude stabilization of an underactuated spacecraft.JournalofAstronautics, 2013,34(1):17～24 (in Chinese))

5莊宇飛,馬廣富,黃海濱. 欠驅(qū)動剛性航天器時間最優(yōu)軌跡規(guī)劃設(shè)計. 控制與決策, 2010,25(10):1469～1473 (Zhuang Y F, Ma G F, Huang H B. Time-optimal motion planning of an underactuated rigid spacecraft.ControlandDecision, 2010,25(10):1469～1473 (in Chinese))

6Benson D A. A gauss pseudospectral transcription for optimal control[PhD Thesis]. Cambridge, Massachusetts：Massachusetts Institute of Technology, 2005

7Tang X J. A Chebyshev-gauss pseudospectral method for solving optimal control problems.ActaAutomaticaSinica, 2015,41(10):1778～1787

8Peterson C. Advances in underactuated spacecraft control[PhD Thesis]. Michigan: University of Michigan, 2016

9Hughes P C. Spacecraft attitude dynamics. New Jersey: Wiley, Hoboken, 1994:17～31

10 Sidi M J. Spacecraft dynamics and control. Cambridge: Cambridge University Press, 1997:25～28

11 Isidori A. Nonlinear control systems II. Berlin:Springer-Verlag, 1999

12 Ge X S, Yi Z G. Optimal control of attitude for coupled-rigid-body spacecraft via Chebyshev-Gauss pseudospectral method.AppliedMathematicsandMechanics, 2017,38(9):1257～1272

13 Berrut J P, Trefethen L N. Barycentric lagrange interpolation.SIAMReview, 2004,46(3):501～517

14 Costa B, Don W S. On the computation of high order pseudo-spectral derivatives.AppliedNumericalMathematics, 2000,33(1):151～159

15 Trefethen L N. Is Gauss quadrature better than Clenshaw-Curtis?SiamReview, 2008,50(1):67～87

16 Waldvogel J. Fast construction of the Fejér and Clenshaw-Curtis quadrature rules.BITNumericalMathematics, 2006,46(1):195～202

17 Krishnan H, Mc clamroch N H, Reyhanoglu M. Attitude stabilization of a rigid spacecraft using two momentum wheel actuators.JournalofGuidanceControl&Dynamics, 1995,18(2):256～263

18 Gong Q, Ross I M, Fahroo F. Costate computation by a chebyshev pseudospectral method.JournalofGuidanceControl&Dynamics, 2010,33(2):623～628