張相宇，張 剛，曹喜濱

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

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

張相宇，張 剛，曹喜濱

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

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

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

0 引 言

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

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

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

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

1 無徑向推力相對運動模型

1.1考慮J2攝動的相對運動模型

(1)

考慮J2攝動時，目標航天器的軌道角速度在LVLH坐標系三個方向的分量為[14]：

(2)

其導數(shù)為：

(3)

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

(4)

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

1.2動力學方程的SDC形式

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

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

(5)

式中：

α

(6)

(7)

(8)

對式(1)中J2攝動加速度項采用一階線性近似：

(9)

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

(10)

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

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

(11)

(12)

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

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

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

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

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

(13)

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

(14)

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

(15)

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

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

(16)

則由式(13)、(14)和(15)構成的最優(yōu)控制問題可通過求解如下的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)

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

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

矩陣Riccati微分方程(17)的解通常得不到解析的表達式，可采用反向數(shù)值積分的方法從tf時刻積分到當前時刻t進行求解，其具體方法可參考下一節(jié)的圖3并忽略其中的M(x,t)項，但該方法在每個控制周期都需要進行積分，具有較大的計算量。因此本節(jié)采用文獻[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)

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

-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控制方法

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

(31)

式中:xf≠0。

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

(32)

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

定義Hamilton函數(shù)：

(33)

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

(34)

得：

(35)

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

(36)

得：

(37)

及其在t=tf時應當滿足的邊界條件為：

(38)

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

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

(39)

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

(40)

對式(39)關于狀態(tài)x求偏導得：

(41)

將式(41)和式(39)代入式(37)，通過對應項相等得：

(42)

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

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

(43)

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

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

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

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

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

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

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

表2 初始時刻的相對位置和相對速度Table 2 The initial relative position and relative velocity

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

(44)

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

圖4和圖6分別為ec=0和ec=0.3時在目標航天器LVLH坐標系下的最優(yōu)轉(zhuǎn)移軌跡在三維坐標系和各個坐標平面內(nèi)的投影。圖5和圖7分別為ec=0和ec=0.3時，作用在追蹤航天器周向和法向的推力加速度曲線。由圖4～7可知，三種控制方法下追蹤航天器均與目標航天器實現(xiàn)了交會。此外將求得的控制量代入原始的非線性方程，當ec=0時，采用解析求解法得到的終端位置誤差為[-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，解析求解法為2.76415×10-3，數(shù)值積分求解法為2.47747×10-3，偽譜法為2.47736×10-3。數(shù)值積分求解法相對偽譜法的最優(yōu)性偏差為0.004%。

當ec=0.3時，采用解析求解法得到的終端位置誤差為[-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，解析求解法為3.76989×10-3，數(shù)值積分求解法為9.72592×10-4，偽譜法為1.03204×10-3。數(shù)值積分求解法相對偽譜法的最優(yōu)性偏差為5.8%。

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

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

4.2非零終端狀態(tài)的編隊重構仿真

本節(jié)針對橢圓軌道下的編隊飛行隊形重構問題進行仿真，取初始時刻目標航天器和追蹤航天器的軌道參數(shù)如表3所示。

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

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

(45)

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

圖8為在主星LVLH坐標系下的最優(yōu)轉(zhuǎn)移軌跡在三維坐標系和各個坐標平面內(nèi)的投影。圖9為作用在追蹤航天器周向和法向的推力加速度。由圖8～9可知，數(shù)值積分求解法與偽譜法的結果非常一致，均可實現(xiàn)隊形的重構任務。在推力加速度的量級上，周向推力小于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，位置誤差仍在分米量級。最優(yōu)性指標方面，數(shù)值積分法為1.87632×10-4，偽譜法為1.80283×10-4，最優(yōu)性指標的偏差為4.1%。此外在仿真過程中，采用數(shù)值積分法的結果作為高斯偽譜法的初值，使得高斯偽譜法能夠快速的收斂。

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

5 結 論

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

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

3) 擴展了SDRE方法使之能夠求解終端狀態(tài)不為零的編隊飛行隊形重構問題，給出了數(shù)值求解策略。該方法不需要初始值猜測，在偏心率達到0.3時相較于偽譜法的最優(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)星編隊控制[J]. 應用科學學報, 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] 雍恩米, 陳磊, 唐國金. 飛行器軌跡優(yōu)化數(shù)值方法綜述[J]. 宇航學報, 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)化的運載器大氣層內(nèi)上升段軌跡快速規(guī)劃方法[J]. 航空學報, 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ù)小推力航天器的深空探測軌道優(yōu)化方法綜述[J]. 力學與實踐, 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]. 宇航學報, 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] 劉利軍, 沈毅, 趙振昊. 基于多項式擬合SDRE的三維導引律設計[J]. 宇航學報, 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] 張銀輝, 楊華波, 江振宇,等. 基于干擾估計的航天器大角度姿態(tài)機動魯棒次優(yōu)控制[J]. 宇航學報, 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

國家自然科學基金(11402062,91438202)；重點實驗室開放基金(HIT.KLOF.MST.201504)

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

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

電話：(0451)86413440- 8507

E-mail: zxyuhit@163.com

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

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

電話：(0451)86413440- 8507

E-mail: zhanggang@hit.edu.cn