石青鑫 戈新生

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

引言

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

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

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

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

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

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

(1)

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

(2)

(3)

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

(4)

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

(5)

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

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

(6)

(7)

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

(8)

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

(9)

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

(10)

(11)

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

(12)

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

(13)

其中：

(14)

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

(15)

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

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

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

(16)

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

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

(17)

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

3 Chebyshev-Gauss偽譜法

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

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

(18)

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

(19)

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

(20)

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

(21)

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

(22)

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

(23)

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

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

(24)

(25)

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

(26)

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

(27)

式中τk為CG點(diǎn).文中采用重心拉格朗日插值,因此微分矩陣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)的時(shí)間區(qū)間為[-1,1),因而狀態(tài)變量的離散表達(dá)式中未包含末端時(shí)間節(jié)點(diǎn)X(τf),對(duì)(16)中的狀態(tài)約束方程在區(qū)間[-1,1]上積分：

(30)

將式(30)左邊展開得：

(31)

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

(32)

(33)

(34)

(35)

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

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

(36)

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

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

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

圖中CGPM代表CG偽譜法.為驗(yàn)證文中姿態(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ù)量級(jí),可以看出該算法求解得到的姿態(tài)運(yùn)動(dòng)以及控制輸入都是合理且有效的.

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

為直觀體現(xiàn)CG偽譜法的優(yōu)勢(shì),文中同時(shí)給出采用Gauss偽譜法計(jì)算的結(jié)果.Gauss偽譜法的配置點(diǎn)為Legendre-Gauss(LG)點(diǎn)；采用Gauss積分來離散性能指標(biāo)函數(shù)中的積分項(xiàng)；應(yīng)用經(jīng)典拉格朗日插值去逼近狀態(tài)和控制變量.Gauss偽譜法同樣采用串行優(yōu)化策略,可行解選取LG點(diǎn)個(gè)數(shù)為6個(gè),最優(yōu)解選取LG點(diǎn)個(gè)數(shù)為60個(gè),狀態(tài)和控制變量的參數(shù)設(shè)置均與CG偽譜法的相同.CG偽譜法和Gauss偽譜法可行解的初始猜測(cè)值均采用MATLAB的隨機(jī)函數(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總的運(yùn)行時(shí)間.表1直觀地表明,CG偽譜法的目標(biāo)函數(shù)的可行解和最優(yōu)解值均比Gauss偽譜法的小,求解時(shí)間比Gauss偽譜法明顯縮短.

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

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

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

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

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

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

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

5 結(jié)論

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

圖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戈新生,陳立群,劉延柱. 帶有兩個(gè)動(dòng)量飛輪剛體航天器的姿態(tài)非完整運(yùn)動(dòng)規(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ū)動(dòng)航天器姿態(tài)穩(wěn)定的分層滑模控制器設(shè)計(jì). 宇航學(xué)報(bào), 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ū)動(dòng)剛性航天器時(shí)間最優(yōu)軌跡規(guī)劃設(shè)計(jì). 控制與決策, 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