張相宇，張 剛，曹喜濱

(哈爾濱工業(yè)大學(xué)衛(wèi)星技術(shù)研究所，哈爾濱150001)

采用SDRE方法的無(wú)徑向推力最優(yōu)軌道控制

張相宇，張 剛，曹喜濱

(哈爾濱工業(yè)大學(xué)衛(wèi)星技術(shù)研究所，哈爾濱150001)

針對(duì)無(wú)徑向推力作用的兩航天器軌道交會(huì)和編隊(duì)衛(wèi)星隊(duì)形重構(gòu)任務(wù)，采用狀態(tài)依賴Riccati方程(SDRE)方法求解了其最優(yōu)軌道控制問(wèn)題。首先考慮J2攝動(dòng)和推力僅存在于追蹤航天器的周向和法向，推導(dǎo)了狀態(tài)依賴配點(diǎn)(SDC)形式的非線性相對(duì)運(yùn)動(dòng)方程。然后針對(duì)終端狀態(tài)為零的軌道交會(huì)問(wèn)題，采用SDRE方法得到了最優(yōu)反饋控制律，并給出了狀態(tài)依賴Riccati微分方程的近似求解策略和數(shù)值求解策略。接著擴(kuò)展了SDRE方法并將其用于終端狀態(tài)不為零的編隊(duì)衛(wèi)星隊(duì)形重構(gòu)問(wèn)題，并給出了相應(yīng)的數(shù)值求解策略。相比于偽譜法等優(yōu)化方法，本文提出的方法不需要初始猜測(cè)值。此外，數(shù)值仿真表明，解析求解Riccati微分方程方法對(duì)于近圓軌道具有較高的精度，數(shù)值計(jì)算方法對(duì)即使偏心率為0.3的橢圓軌道，其最優(yōu)性偏差仍小于6%。

軌道交會(huì)；編隊(duì)重構(gòu)；狀態(tài)依賴Riccati方程(SDRE)；周向推力；軌跡優(yōu)化

0 引 言

為了實(shí)現(xiàn)兩航天器的軌道交會(huì)或編隊(duì)衛(wèi)星的隊(duì)形重構(gòu)任務(wù)，通常采用能夠產(chǎn)生三個(gè)方向推力的推力器配置形式或者采用單個(gè)推力器通過(guò)航天器的姿態(tài)調(diào)整，實(shí)現(xiàn)軌道控制所需的三個(gè)方向的推力。當(dāng)航天器不便于在三個(gè)方向安裝推力器且不便于通過(guò)姿態(tài)調(diào)整(如對(duì)地觀測(cè)衛(wèi)星的對(duì)地軸安裝有相機(jī)，且在軌道機(jī)動(dòng)過(guò)程中仍然保持對(duì)地定向任務(wù))來(lái)實(shí)現(xiàn)三個(gè)方向的推力控制時(shí)，采用無(wú)徑向推力的方式來(lái)進(jìn)行軌道控制是一種較好的選擇。但由于推力方向的特殊性使得系統(tǒng)的控制能力減弱，增強(qiáng)了模型的非線性特性，這給航天器的相對(duì)軌道控制帶來(lái)了新的挑戰(zhàn)。為了解決該問(wèn)題，本文針對(duì)無(wú)徑向推力條件下，航天器軌道交會(huì)或編隊(duì)衛(wèi)星的隊(duì)形重構(gòu)任務(wù)中的最優(yōu)軌道控制問(wèn)題展開(kāi)研究。

對(duì)于僅包含周向推力的平面相對(duì)運(yùn)動(dòng)問(wèn)題，Leonard等[1]提出將兩個(gè)航天器之間的大氣阻力差值等效為周向推力實(shí)現(xiàn)兩個(gè)航天器的編隊(duì)飛行。此后眾多學(xué)者針對(duì)此類周向推力問(wèn)題進(jìn)行了研究[2-3]。對(duì)于其它形式的周向推力，Kumar等[4]設(shè)計(jì)了圓軌道編隊(duì)飛行的線性反饋控制器。針對(duì)三維軌道控制，Starin等[5]設(shè)計(jì)了無(wú)徑向推力下的圓軌道編隊(duì)飛行線性二次型控制器(Linear quadratic regulator，LQR)，Godard等[6]設(shè)計(jì)了非線性控制器進(jìn)行編隊(duì)構(gòu)型保持和隊(duì)形重構(gòu)。Huang等[7-8]研究了兩航天器相對(duì)位置懸停和編隊(duì)重構(gòu)問(wèn)題。以上研究均是基于圓軌道線性相對(duì)運(yùn)動(dòng)方程，少有考慮橢圓軌道的情況，且未考慮J2攝動(dòng)的影響。此外，以上研究中的推力方向均考慮在目標(biāo)航天器LVLH坐標(biāo)系的周向和法向，并不是嚴(yán)格意義上追蹤航天器的周向和法向。

在相對(duì)運(yùn)動(dòng)模型中，目標(biāo)軌道為圓軌道時(shí)常采用C-W方程，橢圓軌道時(shí)常采用T-H方程。在考慮J2攝動(dòng)的相對(duì)運(yùn)動(dòng)模型中，Schweighart和Sedwick[9]采用J2攝動(dòng)的一階近似擴(kuò)展了C-W方程，得到了圓軌道下的線性相對(duì)運(yùn)動(dòng)方程，但未考慮攝動(dòng)中短周期項(xiàng)的影響。Vadali[10]考慮了短周期項(xiàng)的影響得到了平均圓軌道的線性相對(duì)運(yùn)動(dòng)方程。Hamel和Lafontaine[11]采用平均軌道根數(shù)得到了考慮J2攝動(dòng)的橢圓軌道相對(duì)運(yùn)動(dòng)方程。以上均是線性化的相對(duì)運(yùn)動(dòng)模型，對(duì)于非線性模型，Massari等[12]得到了狀態(tài)依賴配點(diǎn)(State-dependent coefficient, SDC)形式的相對(duì)運(yùn)動(dòng)模型，但該模型僅考慮了J2攝動(dòng)對(duì)目標(biāo)航天器的軌道的影響，忽略了兩航天器間J2攝動(dòng)力的差值。Park等[13]雖然考慮了J2攝動(dòng)力的差值，但未考慮J2攝動(dòng)對(duì)目標(biāo)航天器的軌道的影響。Xu和Wang[14]給出了一種精確的相對(duì)運(yùn)動(dòng)方程，但該方程不便于寫成SDC的形式。

本文考慮推力在追蹤航天器的周向和法向，建立SDC形式的非線性相對(duì)運(yùn)動(dòng)方程。擴(kuò)展現(xiàn)有的SDRE方法得到終端狀態(tài)為零和不為零的非線性最優(yōu)反饋控制律，并給出終端狀態(tài)為零的解析計(jì)算方法和兩種情況下的數(shù)值求解策略。分別針對(duì)軌道交會(huì)和編隊(duì)飛行隊(duì)形重構(gòu)問(wèn)題進(jìn)行仿真，以校驗(yàn)本方法對(duì)不同偏心率的主星軌道的求解精度。

1 無(wú)徑向推力相對(duì)運(yùn)動(dòng)模型

1.1考慮J2攝動(dòng)的相對(duì)運(yùn)動(dòng)模型

(1)

考慮J2攝動(dòng)時(shí)，目標(biāo)航天器的軌道角速度在LVLH坐標(biāo)系三個(gè)方向的分量為[14]：

(2)

其導(dǎo)數(shù)為：

(3)

此外，式(1)中ud為作用于追蹤航天器的推力加速度在目標(biāo)航天器LVLH坐標(biāo)系下的表示。當(dāng)考慮推力僅存在于追蹤航天器的周向和法向時(shí)，需對(duì)推力加速度作如下變換：

(4)

式中：uτ和un分別為追蹤航天器的周向和法向推力加速度，Crel為追蹤航天器LVLH坐標(biāo)系到目標(biāo)航天器LVLH坐標(biāo)系的轉(zhuǎn)換矩陣，其具體表達(dá)式將在下一節(jié)進(jìn)行說(shuō)明。

1.2動(dòng)力學(xué)方程的SDC形式

基于SDRE方法的最優(yōu)控制求解法可直接針對(duì)非線性系統(tǒng)進(jìn)行，從而避免線性化帶來(lái)的誤差，但需將系統(tǒng)方程寫成SDC形式，即狀態(tài)變量乘以某個(gè)函數(shù)的形式。式(1)中與角速度相關(guān)的項(xiàng)已是SDC形式，本節(jié)給出二體引力加速度差值項(xiàng)和J2攝動(dòng)加速度差值項(xiàng)的SDC形式以及轉(zhuǎn)換矩陣Crel的二階近似。

對(duì)式(1)中二體引力加速度差值項(xiàng)，采用文獻(xiàn)[23]中的方法可以表示為：

(5)

式中：

α

(6)

(7)

(8)

對(duì)式(1)中J2攝動(dòng)加速度項(xiàng)采用一階線性近似：

(9)

對(duì)于追蹤航天器LVLH坐標(biāo)系到目標(biāo)航天器LVLH坐標(biāo)系的轉(zhuǎn)換矩陣Crel，文獻(xiàn)[24]給出了式(8)所示的二階近似，其中(上標(biāo)“-”)表示歸一化處理后的值，其具體方法可參考文獻(xiàn)[24]。

(10)

其中，系數(shù)A(x)和B(x)見(jiàn)式(11)和式(12)。

需要說(shuō)明的是，本文給出的SDC形式的相對(duì)運(yùn)動(dòng)模型相比文獻(xiàn)[13]進(jìn)一步考慮了J2攝動(dòng)對(duì)主星軌道角速度的影響；相比文獻(xiàn)[12]中的模型，本文更全面的考慮了兩航天器的J2攝動(dòng)差分項(xiàng)。此外本文的模型并沒(méi)有將二階引力差值項(xiàng)進(jìn)行線性化處理。

(11)

(12)

此外，本文的方法也適用于有徑向推力的情況。此時(shí)可在式(4)中加入徑向推力ur項(xiàng)，采用與本文類似的推導(dǎo)方法，仍然可得到形如式(10)的SDC形式的相對(duì)運(yùn)動(dòng)方程，后續(xù)基于SDRE的求解方法的推導(dǎo)過(guò)程也是一致的。

2 終端狀態(tài)為零的SDRE控制方法

在兩航天器的軌道交會(huì)任務(wù)中，要求在終端時(shí)刻兩航天器的相對(duì)位置和相對(duì)速度為零，該問(wèn)題等價(jià)于求解終端時(shí)間給定且終端狀態(tài)為零的最優(yōu)控制問(wèn)題。

2.1狀態(tài)依賴Riccati微分方程的建立

考慮如下的非線性系統(tǒng)：

(13)

其在初始時(shí)刻t0和終端時(shí)刻tf的狀態(tài)滿足：

(14)

以及最優(yōu)控制的性能指標(biāo)為：

(15)

式中：x∈Rn、f(x)∈Rn和u∈Rm為向量，B(x)∈Rn×m為矩陣，Q∈Rn×n、S∈Rn×n和R∈Rm×m分別為對(duì)稱正定矩陣。

將式(13)的系統(tǒng)寫成如下SDC形式：

(16)

則由式(13)、(14)和(15)構(gòu)成的最優(yōu)控制問(wèn)題可通過(guò)求解如下的Riccati微分方程得到[19]：

P(x,t)B(x)R-1BT(x)P(x,t)

(17)

且滿足如下終端條件：

P(x,tf)=S

(18)

則所得的最優(yōu)反饋控制律為：

u(x,t)=-R-1BT(x)P(x,t)x(t)

(19)

此外，由于該方法為反饋控制，因此對(duì)初始偏差也有一定的適應(yīng)能力。

2.2狀態(tài)依賴Riccati微分方程的近似求解策略

矩陣Riccati微分方程(17)的解通常得不到解析的表達(dá)式，可采用反向數(shù)值積分的方法從tf時(shí)刻積分到當(dāng)前時(shí)刻t進(jìn)行求解，其具體方法可參考下一節(jié)的圖3并忽略其中的M(x,t)項(xiàng)，但該方法在每個(gè)控制周期都需要進(jìn)行積分，具有較大的計(jì)算量。因此本節(jié)采用文獻(xiàn)[19]提出的一種近似求解策略。

首先求解如下的代數(shù)Riccati方程：

Pss(x)A(x)+ATPss(x)-

Pss(x)B(x)R-1(x)Pss(x)+Q=0

(20)

將式(20)代入式(17)，可得：

AT(x)(P(x,t)-Pss(x))-

P(x,t)B(x)R-1BT(x)P(x,t)+

Pss(x)B(x)R-1BT(x)Pss(x)

(21)

與文獻(xiàn)[19]所不同的是，為了避免在后續(xù)對(duì)P(x,t)-Pss(x)求逆的過(guò)程中出現(xiàn)奇異，本文并不直接求取式(20)的正定解，而采用文獻(xiàn)[25]的方法求取其負(fù)定解。通過(guò)求取如下方程

-Pn(x)A(x)-ATPn(x)-

Pn(x)B(x)R-1(x)Pn(x)+Q=0

(22)

(23)

(24)

式(21)可以表示為：

B(x)R-1BT(x)

(25)

其終端條件為：

(26)

在方程(25)中令

(27)

則式(25)的解為：

(28)

式中:E為如下代數(shù)Lyapunov方程的解

(29)

至此，求得Riccati微分方程(17)的解為：

(30)

綜上，解析求解Riccati微分方程(17)的步驟如圖2所示。

3 終端狀態(tài)不為零的SDRE控制方法

在編隊(duì)飛行的隊(duì)形重構(gòu)任務(wù)中，繞飛航天器在給定的時(shí)間內(nèi)相對(duì)主航天器從一個(gè)相對(duì)位置到達(dá)另一個(gè)相對(duì)位置從而實(shí)現(xiàn)新的構(gòu)型。此時(shí)為了滿足新的繞飛條件，終端狀態(tài)通常不為零。假設(shè)追蹤航天器的相對(duì)狀態(tài)在初始時(shí)刻t0和終端時(shí)刻tf分別為：

(31)

式中:xf≠0。

取最優(yōu)性能指標(biāo)：

(32)

與式(15)不同的是，該性能指標(biāo)中與狀態(tài)x相關(guān)的項(xiàng)不再是標(biāo)準(zhǔn)的二次型的形式，因此第2節(jié)關(guān)于終端狀態(tài)為零的SDRE方法并不能直接用于本節(jié)問(wèn)題的求解。本節(jié)基于動(dòng)態(tài)規(guī)劃方法中的HJB方程，將文獻(xiàn)[19]中標(biāo)準(zhǔn)的SDRE方法擴(kuò)展到終端狀態(tài)不為零的情形。

定義Hamilton函數(shù)：

(33)

其中,J*(x,t)稱為最優(yōu)值函數(shù)，由極大值原理：

(34)

得：

(35)

將式(35)代入式(36)，再由連續(xù)系統(tǒng)的HJB方程：

(36)

得：

(37)

及其在t=tf時(shí)應(yīng)當(dāng)滿足的邊界條件為：

(38)

由于邊界條件包含x的二次、一次、零次多項(xiàng)式的形式，因此假設(shè)最優(yōu)值函數(shù)J*(x,t)具有如下與邊界條件一致的形式：

xT(t)M(x,t)+N(x,t)

(39)

為保證最優(yōu)值函數(shù)滿足邊界條件(38)，其中P、M、N在終端時(shí)刻應(yīng)當(dāng)滿足：

(40)

對(duì)式(39)關(guān)于狀態(tài)x求偏導(dǎo)得：

(41)

將式(41)和式(39)代入式(37)，通過(guò)對(duì)應(yīng)項(xiàng)相等得：

(42)

其中，P、M、N在終端時(shí)刻滿足式(40)的邊界條件，相應(yīng)的反饋控制可以表示為：

u=-R-1(x)BT(x)(P(x,t)x(t)+M(x,t))

(43)

由式(40)和式(42)可知：當(dāng)xf=0時(shí)M(x,t)≡0，此時(shí)所得的控制律(43)與第2節(jié)的控制律(19)是一致的，因此第2節(jié)終端狀態(tài)為零的SDRE方法是本節(jié)所得結(jié)果的特殊情形。此外，在式(42)中M(x,t)是P(x,t)的函數(shù)，因此很難求得M(x,t)的解析表達(dá)式。圖3給出一種數(shù)值求解策略：在每一個(gè)正向積分步驟k中，先采用逆向數(shù)值積分的方法計(jì)算出從tf時(shí)刻到當(dāng)前時(shí)刻tk的Riccati矩陣P(xk,tk)和M(xk,tk)，然后通過(guò)式(43)計(jì)算出本步長(zhǎng)內(nèi)的控制量uk。在此過(guò)程中目標(biāo)航天器的軌道參數(shù)可通過(guò)事先求得的離散軌道參數(shù)插值得到。此外，在逆向積分過(guò)程中狀態(tài)x在tf到tk間的值并不知曉，因此在逆向積分過(guò)程中將兩航天器的相對(duì)狀態(tài)固定為tk時(shí)的值x(tk)。此外，由于N(x,t)與M(x,t)和P(x,t)不耦合且并沒(méi)有出現(xiàn)在式(42)的控制中，因此在逆向數(shù)值積分式(42)時(shí)，僅積分前兩式以減小計(jì)算量。雖然該方法在每個(gè)控制周期都需要進(jìn)行數(shù)值積分運(yùn)算具有較大的計(jì)算量，但避免了優(yōu)化過(guò)程中的初值猜測(cè)，因此仍然具有較大的實(shí)用價(jià)值。

4 數(shù)值仿真及結(jié)果分析

本節(jié)分別針對(duì)近距離軌道交會(huì)和編隊(duì)飛行的隊(duì)形重構(gòu)問(wèn)題進(jìn)行仿真校驗(yàn)。其中軌道交會(huì)問(wèn)題對(duì)應(yīng)終端狀態(tài)為零的情況，重構(gòu)問(wèn)題對(duì)應(yīng)終端狀態(tài)不為零的情況。在交會(huì)問(wèn)題中將比較解析求解法、數(shù)值積分求解法和高斯偽譜法三種方法的優(yōu)化結(jié)果，在重構(gòu)問(wèn)題中將比較數(shù)值積分求解法和高斯偽譜法的優(yōu)化結(jié)果。高斯偽譜法采用佛羅里達(dá)大學(xué)開(kāi)發(fā)的GPOPS優(yōu)化軟件。仿真中用到的基本物理參數(shù)為：地球半徑Re=6378.13 km，地球引力常數(shù)為μ=3.986004×105km3/s2，攝動(dòng)系數(shù)J2=1.082629×10-3。

4.1終端狀態(tài)為零的軌道交會(huì)仿真

為了比較本文提出的算法對(duì)不同偏心率的目標(biāo)航天器軌道的有效性，分別針對(duì)偏心率ec=0和ec=0.3兩種情況進(jìn)行了仿真。初始時(shí)刻目標(biāo)航天器和追蹤航天器的軌道參數(shù)如表1所示。

表1 初始時(shí)刻兩航天器的軌道參數(shù)Table 1 The initial orbit elements of two spacecrafts

表1中,Hc為目標(biāo)航天器的近地點(diǎn)高度，下標(biāo)“c”表示目標(biāo)航天器的軌道參數(shù)，下標(biāo)“d”表示追蹤航天器的軌道參數(shù)。設(shè)置軌道交會(huì)的初始時(shí)刻t0=0，終端時(shí)刻tf=1.3Tc，其中Tc為目標(biāo)航天器的初始軌道周期。根據(jù)以上軌道參數(shù)可得初始的相對(duì)位置和相對(duì)速度在LVLH坐標(biāo)系下的值如表2所示。

表2 初始時(shí)刻的相對(duì)位置和相對(duì)速度Table 2 The initial relative position and relative velocity

此外，對(duì)于燃料最優(yōu)問(wèn)題，性能指標(biāo)函數(shù)(15)中應(yīng)取Q=0且S為無(wú)窮大，但在解析求解法中為了保證Riccati方程解的存在性，并且考慮到位置和速度在數(shù)量級(jí)上的差別，取

(44)

在數(shù)值積分方法中不需要對(duì)Q進(jìn)行限制，因此取Q=0，但S仍與式(44)一致。在偽譜法中由于方法原理的不同，可以直接在性能指標(biāo)中不包含Q和S。

圖4和圖6分別為ec=0和ec=0.3時(shí)在目標(biāo)航天器LVLH坐標(biāo)系下的最優(yōu)轉(zhuǎn)移軌跡在三維坐標(biāo)系和各個(gè)坐標(biāo)平面內(nèi)的投影。圖5和圖7分別為ec=0和ec=0.3時(shí)，作用在追蹤航天器周向和法向的推力加速度曲線。由圖4～7可知，三種控制方法下追蹤航天器均與目標(biāo)航天器實(shí)現(xiàn)了交會(huì)。此外將求得的控制量代入原始的非線性方程，當(dāng)ec=0時(shí)，采用解析求解法得到的終端位置誤差為[-0.03580,0.34287,-0.02413]Tm，終端速度誤差為[-0.24981,6.19937,7.68697]T×10-5m/s；采用數(shù)值積分求解法得到的終端位置誤差為[-0.03827,0.35629,-0.02517]Tm，終端速度誤差為[-0.95065,5.49735,7.70494]T×10-5m/s，位置誤差在分米量級(jí)。對(duì)于性能指標(biāo)J，解析求解法為2.76415×10-3，數(shù)值積分求解法為2.47747×10-3，偽譜法為2.47736×10-3。數(shù)值積分求解法相對(duì)偽譜法的最優(yōu)性偏差為0.004%。

當(dāng)ec=0.3時(shí)，采用解析求解法得到的終端位置誤差為[-0.05341,0.12008,-0.00334]Tm，終端速度誤差為[-4.96389,1.53857,0.03782]T×10-3m/s；采用數(shù)值積分求解法得到的終端位置誤差為[-0.11342,0.18104,-0.01197]Tm，終端速度誤差為[-3.86410,5.62476,1.41895]T×10-5m/s。位置誤差仍在分米量級(jí)。對(duì)于性能指標(biāo)J，解析求解法為3.76989×10-3，數(shù)值積分求解法為9.72592×10-4，偽譜法為1.03204×10-3。數(shù)值積分求解法相對(duì)偽譜法的最優(yōu)性偏差為5.8%。

由以上結(jié)果可知，隨著偏心率ec的增大，解析計(jì)算方法和數(shù)值積分法的結(jié)果與偽譜法的偏差逐漸增大。數(shù)值積分法的偏差主要來(lái)自于每次逆向積分計(jì)算Riccati矩陣的過(guò)程中，將系數(shù)矩陣A(x,t),B(x,t)中的狀態(tài)量x視為與當(dāng)前時(shí)刻相等的常值，但實(shí)際中隨著時(shí)間逐漸接近終端時(shí)刻，狀態(tài)x是逐漸接近于零的。解析計(jì)算方法的偏差大于數(shù)值積分方法的偏差，主要原因是在解析計(jì)算Riccati矩陣的過(guò)程中，不僅將系數(shù)矩陣A(x,t),B(x,t)中的狀態(tài)量x視為與當(dāng)前時(shí)刻相等的常值，且將目標(biāo)航天器的軌道參數(shù)也視為與當(dāng)前時(shí)刻相等的常數(shù)。由式(11)和式(12)可知，對(duì)系數(shù)A21和A22的影響包含兩航天器的相對(duì)位置x和目標(biāo)航天器的軌道參數(shù)兩部分，其中與目標(biāo)航天器軌道角速度ωc相關(guān)的項(xiàng)是周期性變化的，且幅值隨著偏心率ec的增大而增大。因此，在解析法求解過(guò)程的每一步將A(x,t),B(x,t)視為常數(shù)，其所得結(jié)果的最優(yōu)性必然會(huì)隨著偏心率ec的增大而降低。

雖然解析求解法的最優(yōu)性會(huì)隨著偏心率ec的增大而降低，但仿真過(guò)程中將其結(jié)果作為偽譜法的初值，使得偽譜法的收斂速度大大提高。

4.2非零終端狀態(tài)的編隊(duì)重構(gòu)仿真

本節(jié)針對(duì)橢圓軌道下的編隊(duì)飛行隊(duì)形重構(gòu)問(wèn)題進(jìn)行仿真，取初始時(shí)刻目標(biāo)航天器和追蹤航天器的軌道參數(shù)如表3所示。

表3 初末時(shí)刻兩航天器的軌道參數(shù)(非零終端)Table 3 Initial and final orbit elements of the two spacecraft(non-zero terminal condition)

仿真初始時(shí)刻t0=0，終端時(shí)刻tf=1.2Tc。需要說(shuō)明的是，上表中下標(biāo)“c0”和“ct”分別表示主航天器在初始和終端時(shí)刻的軌道參數(shù)。根據(jù)以上參數(shù)，可以算出LVLH坐標(biāo)系下的初末相對(duì)位置和相對(duì)速度如下：

(45)

由于在J2攝動(dòng)下并不存在嚴(yán)格的周期繞飛條件，但通過(guò)以上方法得到的相對(duì)位置和相對(duì)速度在一個(gè)主航天器軌道周期內(nèi)的偏差較小可近似認(rèn)為是周期軌道。Q、S和R的取值與上一節(jié)一致。

圖8為在主星LVLH坐標(biāo)系下的最優(yōu)轉(zhuǎn)移軌跡在三維坐標(biāo)系和各個(gè)坐標(biāo)平面內(nèi)的投影。圖9為作用在追蹤航天器周向和法向的推力加速度。由圖8～9可知，數(shù)值積分求解法與偽譜法的結(jié)果非常一致，均可實(shí)現(xiàn)隊(duì)形的重構(gòu)任務(wù)。在推力加速度的量級(jí)上，周向推力小于4×10-4m/s2，法向推力小于1×10-4m/s2。此外將求得的控制量代入原始的非線性方程，數(shù)值積分求解法得到的終端位置誤差為[-0.03490,-0.29369,-0.00802]Tm，終端速度誤差為[-3.04687,-1.41242,-0.77204]T×10-5m/s，位置誤差仍在分米量級(jí)。最優(yōu)性指標(biāo)方面，數(shù)值積分法為1.87632×10-4，偽譜法為1.80283×10-4，最優(yōu)性指標(biāo)的偏差為4.1%。此外在仿真過(guò)程中，采用數(shù)值積分法的結(jié)果作為高斯偽譜法的初值，使得高斯偽譜法能夠快速的收斂。

雖然本文的方法需要進(jìn)行數(shù)值積分，但相較于偽譜法，本文提出的方法并不需要進(jìn)行初始猜測(cè)，因此具有很好的實(shí)用性。

5 結(jié) 論

1) 針對(duì)推力在追蹤航天器的周向和法向的情況，建立了SDC形式的考慮J2攝動(dòng)的非線性相對(duì)運(yùn)動(dòng)方程。

2) 基于SDRE方法給出了終端狀態(tài)為零的軌道交會(huì)問(wèn)題的近似解析求解方法。在近圓軌道下具有較好的最優(yōu)性，盡管隨著目標(biāo)航天器偏心率的增大最優(yōu)性降低，但仍可作為偽譜法等直接優(yōu)化算法的初值。

3) 擴(kuò)展了SDRE方法使之能夠求解終端狀態(tài)不為零的編隊(duì)飛行隊(duì)形重構(gòu)問(wèn)題，給出了數(shù)值求解策略。該方法不需要初始值猜測(cè)，在偏心率達(dá)到0.3時(shí)相較于偽譜法的最優(yōu)性偏差仍小于6%。

[1] Leonard C L, Hollister W M, Bergmann E V. Orbital formationkeeping with differential drag[J]. Journal of Guidance, Control, and Dynamics, 1989, 12(1): 108-113.

[2] Bevilacqua R, Romano M. Rendezvous maneuvers of multiple spacecraft using differential drag underJ2perturbation[J]. Journal of Guidance, Control, and Dynamics, 2008, 31(6): 1595-1607.

[3] 陳志明, 王惠南, 劉海穎. 基于大氣阻力的微小衛(wèi)星編隊(duì)控制[J]. 應(yīng)用科學(xué)學(xué)報(bào), 2010, 28(2): 209-215. [Chen Zhi-ming, Wang Hui-nan, Liu Hai-ying. Satellite formation fly control based on atmospheric drag[J]. Journal of Applied Sciences, 2010, 28(2): 209-215.]

[4] Kumar K D, Bang H C, Tahk M J. Satellite formation flying using along-track thrust[J]. Acta Astronautica, 2007, 61(7-8): 553-564.

[5] Starin S R, Yedavalli R K, Sparks A G. Spacecraft formation flying maneuvers using linear quadratic regulation with no radial axis inputs[C]. AIAA Guidance, Navigation, and Control Conference and Exhibit, Montreal, Canada, August 6-9, 2001.

[6] Godard, Kumar K D, Zou A. Robust stationkeeping and reconfiguration of underactuated spacecraft formations[J]. Acta Astronautica, 2014, 105(2): 495-510.

[7] Huang X, Yan Y, Zhou Y. Nonlinear control of underactuated spacecraft hovering[J]. Journal of Guidance, Control, and Dynamics, 2015, 39(3): 685-694.

[8] Huang X, Yan Y, Zhou Y. Analytical solutions to optimal underactuated spacecraft formation reconfiguration[J]. Advances in Space Research, 2015, 56(10): 2151-2166.

[9] Schweighart S A, Sedwick R J. High-fidelity linearizedJ2model for satellite formation flight[J]. Journal of Guidance Control and Dynamics, 2002, 25(6): 1073-1080.

[10] Vadali S R. Model for linearized satellite relative motion about aJ2-perturbed mean circular orbit[J]. Journal of Guidance Control and Dynamics, 2009, 32(5): 1687-1691.

[11] Hamel J, Lafontaine J D. Linearized dynamics of formation flying spacecraft on aJ2-perturbed elliptical orbit[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(6): 1649-1658.

[12] Massari M, Bernelli-Zazzera F, Canavesi S. Nonlinear control of formation flying with state constraints[J]. Journal of Guidance Control and Dynamics, 2012, 35(6): 1919-1925.

[13] Park H, Park S, Choi K. Satellite formation reconfiguration and station-keeping using state-dependent Riccati equation technique[J]. Aerospace Science and Technology, 2011, 15(6): 440-452.

[14] Xu G, Wang D. Nonlinear dynamic equations of satellite relative motion around an oblate earth[J]. Journal of Guidance Control and Dynamics, 2008, 31(5): 1521-1524.

[15] 雍恩米, 陳磊, 唐國(guó)金. 飛行器軌跡優(yōu)化數(shù)值方法綜述[J]. 宇航學(xué)報(bào), 2008, 29(2): 397-406. [Yong En-mi, Chen Lei, Tang Guo-jin. A survey of numerical methods for trajectory optimization of spacecraft[J]. Journal of Astronautics, 2008, 29(2): 397-406.]

[16] 崔乃剛, 黃盤興, 路菲, 等. 基于混合優(yōu)化的運(yùn)載器大氣層內(nèi)上升段軌跡快速規(guī)劃方法[J]. 航空學(xué)報(bào), 2015, 36(6): 1915-1923. [Cui Nai-gang, Huang Pan-xing, Lu Fei, et al. A hybrid optimization approach for rapid endo-atmospheric ascent trajectory planning of launch vehicles[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(6): 1915-1923.]

[17] 李俊峰, 蔣方華. 連續(xù)小推力航天器的深空探測(cè)軌道優(yōu)化方法綜述[J]. 力學(xué)與實(shí)踐, 2011, 33(3): 1-6. [Li Jun-feng, Jiang Fang-hua. Survey of low-thrust trajectory optimization methods for deep space exploration[J]. Mechanics in Engineering, 2011, 33(3): 1-6.]

[18] ?imen T. Survey of state-dependent Riccati equation in nonlinear optimal feedback control synthesis[J]. Journal of Guidance, Control, and Dynamics, 2012, 35(4): 1025-1047.

[19] Heydari A, Balakrishnan S N. Path planning using a novel finite horizon suboptimal controller[J]. Journal of Guidance Control and Dynamics, 2013, 36(4): 1210-1214.

[20] 張軍, 徐世杰. 基于SDRE方法的撓性航天器姿態(tài)控制[J]. 宇航學(xué)報(bào), 2008, 29(1): 138-144. [Zhang Jun, Xu Shi-jie. Control of flexible spacecraft via state-dependent Riccati equation technique[J]. Journal of Astronautics, 2008, 29(1): 138-144.]

[21] 劉利軍, 沈毅, 趙振昊. 基于多項(xiàng)式擬合SDRE的三維導(dǎo)引律設(shè)計(jì)[J]. 宇航學(xué)報(bào), 2010, 31(1): 87-92. [Liu Li-jun, Shen Yi, Zhao Zhen-hao. Three-dimensional missile guidance law design based on polynomial fitting of SDRE [J]. Journal of Astronautics, 2010, 31(1): 87-92.]

[22] 張銀輝, 楊華波, 江振宇,等. 基于干擾估計(jì)的航天器大角度姿態(tài)機(jī)動(dòng)魯棒次優(yōu)控制[J]. 宇航學(xué)報(bào), 2015, 36(10): 1148-1154. [Zhang Yi-hui, Yang Hua-bo, Jiang Zhen-yu, et al. Spacecraft large angle attitude maneuver robust suboptimal control based on disturbance estimation [J]. Journal of Astronautics, 2015, 36(10): 1148-1154.]

[23] Franzini G, Innocenti M. Nonlinear H-infinity control of relative motion in space via the state-dependent Riccati equations[C]. IEEE 54th Conference on Decision and Control (CDC), Osaka, Japan,December 15-18, 2015.

[24] Sengupta P. Elliptic rendezvous in the chaser satellite frame[J]. The Journal of the Astronautical Sciences, 2012, 59(1-2): 216-236.

[25] Nguyen T, Gajic Z. Solving the matrix differential Riccati equation: a Lyapunov equation approach[J]. IEEE Transactions on Automatic Control, 2010, 55(1): 191-194.

OptimalOrbitalControlwithoutRadialThrustBasedonSDREMethod

ZHANG Xiang-yu, ZHANG Gang, CAO Xi-bin

(Research Center of Spacecraft Technology, Harbin Institute of Technology, Harbin 150001, China)

The optimal control problem of spacecraft rendezvous and formation reconfiguration without radial thrust is investigated by using the state dependent Riccati equation (SDRE) method. Considering the thrusters installed in the circumferential and normal directions, the nonlinear relative motion equations influenced byJ2perturbation are derived and presented via the state-dependent coefficient (SDC). Then the optimal feedback control law is derived by using the SDRE method, and an approximate analytical solution of the state dependent Riccati differential equation is obtained. Moreover, a numerical version of the obtained analytical solution is extended to solve the formation reconfiguration problem. Compared with the Gauss pseudospectral spectral method, the proposed method does not need an initial guess. Numerical results show that the analytical solution of the Riccati differential equation method has high precision for the near circular orbits, and the error of the optimal index of the numerical method is less than 6% even for the elliptical orbits with an eccentricity approaching to 0.3.

Orbital rendezvous; Formation reconfiguration; State dependent Riccati equation (SDRE); Circumferential thrust; Trajectory optimization

V448.234

A

1000-1328(2017)10- 1057- 11

10.3873/j.issn.1000-1328.2017.10.006

2016- 09- 12;

2017- 04- 21

國(guó)家自然科學(xué)基金(11402062,91438202)；重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金(HIT.KLOF.MST.201504)

張相宇(1986-)，男，博士生，主要從事航天器導(dǎo)航、制導(dǎo)與控制研究。

通信地址：黑龍江省哈爾濱市南崗區(qū)一匡街2號(hào)哈爾濱工業(yè)大學(xué)科學(xué)園B3棟507室(150080)

電話：(0451)86413440- 8507

E-mail: zxyuhit@163.com

張剛(1983-)，男，博士，副教授，主要從事航天器導(dǎo)航、制導(dǎo)與控制研究。本文通信作者。

通信地址：黑龍江省哈爾濱市南崗區(qū)一匡街2號(hào)哈爾濱工業(yè)大學(xué)科學(xué)園B3棟507室(150080)

電話：(0451)86413440- 8507

E-mail: zhanggang@hit.edu.cn