一種基于多終端約束的最優制導方法

李超兵，王晉麟，李海

北京航天自動控制研究所，北京 100854

?

一種基于多終端約束的最優制導方法

李超兵*，王晉麟，李海

北京航天自動控制研究所，北京 100854

在航天器主發動機推力大小不可調的前提下，針對5個終端約束下傳統迭代制導小角度修正假設的不足，對一種基于多終端約束的最優制導方法進行了研究。在入軌點軌道坐標系下建立航天器的最優控制模型，對橫截條件方程組直接進行迭代求解獲得制導角度指令，在此基礎上，通過對開關機點進行優化以減小未被滿足的終端位置約束的影響；進一步，推導了地心慣性系下等效的5個終端約束，并通過引入權重因子來提高制導方程數值求解的精度。標準條件下的仿真結果表明，所提制導方法與傳統迭代制導相比，未被滿足的終端位置約束精度提高了159.535 5m，而其余5個終端約束幾乎不受影響；蒙特卡羅打靶仿真結果表明，所提制導方法對航天器初始位置和速度偏差具有一定的適用性。

航天器變軌；迭代制導；最優控制；坐標轉換；開關機點優化；權重調整

航天器在空間中為完成預定的任務，需要通過制導來完成不同軌道之間的轉移。由于航天器的發動機推力大小有限，因此通常需要選取合適的發動機開關機點以完成制導[1-4]。

傳統的制導方法包括攝動制導、迭代制導等，由于攝動制導存在射前裝訂數據復雜，入軌精度差等缺點，為進一步提高制導精度，以迭代制導代表的自適應制導方案逐漸為許多航天器所采納[5]。

迭代制導最初應用于美國的土星五號運載器上[6]，國內韓祝齋[7]在最優控制原理的框架下對迭代制導的方程進行了描述；陳新民和余夢倫[8]推導了多級迭代制導方程，茹家欣[9]給出了一種更為簡潔的迭代制導方法。

傳統的迭代制導針對航天器主發動機推力大小不可調的情況，通常考慮5個終端約束，包括兩個方向的位置約束和3個方向的速度約束，首先求解滿足速度約束的控制角，然后在位置約束引起的角度變化為小量的假設下對控制角進行修正[6-9]，但在某些變軌情形下，這一假設不再成立，此時傳統迭代制導無法適用。

現有文獻對迭代制導的改進包括控制對象模型復雜性和利用最優控制原理推導的方程形式等方面。文獻[10]考慮了受復雜約束情形下的火箭多級發射上升段制導；文獻[11-12]考慮了大氣環境對飛行器制導的影響，在此基礎上利用最優控制原理對制導方程進行推導；文獻[13-14]從控制角的形式等不同角度對迭代制導進行了改進，但均未擺脫傳統迭代制導小角度修正假設的限制；文獻[15]直接考慮了軌道要素約束下的制導方法，但方程形式不夠直觀，特別是對乘子變量的消去處理形式復雜，且數值求解時對終端約束方程的權重調整不變。此外，文獻[9]雖然直接以推力方向的3分量為控制變量進行了推導，但同樣認為位置約束引起的修正為小量。

考慮到上述文獻的不足，本文提出了一種基于多終端約束的最優制導方法。結合制導特點對發動機關機點進行優化，并引入權重因子，通過直接迭代求解橫截條件方程來提高制導精度，并運用數值仿真驗證了該方法的有效性。

1　最優制導方法研究

1.1問題描述

設航天器在某一軌道完成飛行任務后，進入一條自由滑行的轉移軌道L1，此后需要選擇合適的開關機點進行主發動機的開關機制導，從而進入到新的目標軌道L2，且制導結束后航天器需要以指定的精度在L2軌道上自由滑行至預定的下一次任務點Dfinal。航天器主發動機推力大小不可調。

上述問題涉及到兩個方面的需求，一是需要選取合適的開關機點，二是在開關機點確定后，設計合適的制導算法。

1.2開關機點優化

航天器利用主發動機變軌本質上是對脈沖變軌的近似，設脈沖變軌的時刻為tc，由于主發動機的推力和比沖恒定，因此通常的做法是在L1和L2軌道的“交點”處，利用齊奧爾科夫斯基公式計算點火時間Δt，則主發動機的開機時刻為tc-Δt/2，關機時刻為tc+Δt/2。

上述方法所選取的開關機點只是一個初選方案，為保證最終航天器到達Dfinal的精度，該精度直接取決于航天器制導的精度，因此還需要對開關機點進行優化。

通常制導的坐標系有慣性坐標系(包括地心慣性系和發射慣性系)和入軌點軌道坐標系兩類，考慮到入軌點軌道坐標系在幾何方面更能直觀地分析航天器的運動，也更能直觀地反映航天器在制導結束后終端約束的滿足情況，而航天器制導的精度正是通過終端約束來體現的，因此本文主要對入軌點軌道坐標系下制導的情形進行分析。

在入軌點軌道坐標系(ocf坐標系)下，終端約束通常選為Yocf和Zocf方向的位置約束，以及Xocf、Yocf和Zocf三個方向的速度約束，則目標點和航天器制導結束后到達的實際點如圖1所示。

圖1　入軌點軌道坐標系下的目標點和實際點Fig.1　Object point and real position in the orbital coordinate system of orbit-insertion point

由圖1可知，由于未對Xocf方向的位置進行約束，制導結束后Xocf方向存在位置偏差ΔXocf，如果將目標點向“靠近”實際點的方向“滑行”一些(即目標關機時間比初始計算值多數一些)，那么最終的ΔXocf將會減小一些。具體來說，若ΔXocf為正，則目標點正向“滑行”一些；若ΔXocf為負，則目標點反向“滑行”一些。

與開機點(起始點)相比，關機點(目標點)和誤差之間的關系要更為明顯。為便于分析，對初選的開機點一般不作調整(即認為初選的開機點足夠準確)，而將目標點正向或者反向“滑行”一些，即將L2軌道從“交點”處開始相應的“滑行”時間在Δt/2的基礎上，適當增加或者減小一些，因此，這實際上相當于引入了一個新的控制變量，即“滑行”時間的調整量。

由于一旦對目標點進行調整，則整個入軌點軌道坐標系也會發生變化，而制導指令和最后的誤差都是在入軌點軌道坐標系中計算的，相應的制導流程也會發生變化，因此調整前后的ΔXocf之間存在一個復雜的關系，無法通過簡單的數學推導來得到具體的調整量，需要進行簡化分析。

考慮到在入軌點軌道坐標系中，目標點的速度矢量位于OXocfYocf平面內，一個自然的想法是“滑行”時間的調整量與位置偏差ΔXocf和Xocf方向的速度分量Vxocff存在關系。由于ΔXocf相對于L2軌道的高度來說是小量(迭代制導通常假設其余終端約束滿足的前提下，Xocf方向的位置約束會盡量滿足)，在目標點“靠近”實際點的過程中速度矢量變化很小，且由前面的分析可知，“滑行”時間的調整量應該與ΔXocf正相關，因此不妨將ΔXocf/Vxocff作為“滑行”時間調整的一個基準偏差，根據該偏差來反饋調節“滑行”時間。設“滑行”時間為tcoast，tcoast的調整量正比于ΔXocf/Vxocff，并且可以在前面加一個反饋系數k>0，即“滑行”時間的調整量Δtcoast=kΔXocf/Vxocff，當把ΔXocf/Vxocff看作是誤差時，該表達式類似于經典控制理論中的比例反饋修正，根據前面的分析，開機點(起始點)不變，“滑行”時間的調整如圖2所示。

圖2　“滑行”時間調整示意Fig.2　Sketch map of coasting time adjustment

由圖2，具體的“滑行”時間調整步驟如下：

1)利用齊奧爾科夫斯基公式獲得初選的起始點和目標點，即L2軌道從“交點”處正向的初始“滑行”時間tcoast=Δt/2；

2)利用制導算法對設定的起始點和目標點進行制導計算，獲得最后的位置偏差ΔXocf和目標點x方向的速度分量Vxocff；

3)計算“滑行”時間的調整量Δtcoast=kΔXocf/Vxocff；

此外，由上述計算步驟可知，在每一次迭代計算過程中，都需要進行一次制導算法的運算來獲得最后的位置偏差，制導算法應選擇在入軌點軌道坐標系下進行解算。

1.3最優制導算法

傳統的迭代制導從最優控制原理出發，但并沒有對獲得的方程組直接求解，而是作了一定的簡化假設，特別是假設了指令角的形式，在此基礎上推導出指令角的相關參數。本文在迭代制導的基礎上，介紹一種直接基于最優控制原理的制導方案，通過直接考慮終端約束來獲得相關方程，對相關方程進行求解來獲得姿態指令角，整個制導過程的思想更為清晰和明確，且同樣滿足最終的入軌精度需求。

經過上述處理后，可得到無量綱化的航天器運動方程為

(1)

式中：r為航天器位置矢量；v為航天器速度矢量；g為地球引力加速度矢量；T為航天器的推力大小；ms為恒定的秒流量；m0為初始質量；u為推力方向矢量。

下面建立航天器制導的最優控制模型，系統方程即為式(1)的狀態空間形式，取狀態變量為

(2)

控制量即推力方向u，即

(3)

則狀態方程為

(4)

性能指標為燃料最省，亦即推力飛行時間最短，即

(5)

約束條件為推力方向的單位化約束，即

(6)

對于邊界條件，初始時刻的狀態給定，即

(7)

E1=Xf(2)-Yocff=0，E2=Xf(3)-Zocff=0，

(8)

列寫哈密頓函數：

(9)

式中：λr、λv和λu為協態變量，則相應的伴隨方程為

(10)

由式(10)可解得

(11)

式中：λr0和λv0分別為λr和λv的初值。

控制方程為

(12)

考慮約束式(6)，可解得

(13)

于是由系統方程解得最優終端狀態為

(14)

式(14)中的積分運算可利用數值插值積分公式進行。

終端約束對應的橫截條件為

(15)

(16)

(17)

由于：

(18)

(19)

(20)

(21)

式(8)、式(17)和式(21)組成了所有7個橫截條件方程，結合協變量和最優狀態量的表達式(10)和式(14)，則橫截條件方程組含λ0和tf7個變量，由于最優終端狀態可表示為λ0和tf的函數，而最優終端狀態需要滿足式(8)中的終端約束，因此可通過數值迭代求解出λ0和tf。在解得λ0和tf后，當前時刻的控制量即為

(22)

需要說明的是，上述方法通過直接迭代求解橫截條件方程組來獲得協態變量初值，而沒有像傳統迭代制導方法那樣對指令角度做簡化假設，因此理論上來說可以獲得最優解。在具體應用時，數值求解的精度決定了獲得控制解的好壞，這一點主要體現在終端約束滿足的效果上。下面對此進一步進行分析。

2　終端約束的權重調整

(23)

(24)

(25)

(26)

(27)

由于入軌點軌道坐標系本身的特點，Yocff是一相當大的數，而理論上Zocff=0，因此直接以式(23)～(27)作為終端約束進行控制量的求解時，由于E1與其余4個約束的量級相差較大，直接迭代求解時會使得終端約束E1實際的值較大，因此需要對原始的終端約束條件進行轉化。

考慮到通常情況下在地心慣性坐標系下目標點3個方向的位置分量大致處于同一量級，因此，如果能夠將入軌點軌道坐標系下的位置約束等效轉化為地心慣性坐標系下的位置約束，那么理論上來說對于迭代求解協態變量的初值沒有任何影響，但從數值角度來說將更為容易處理。

(28)

則地心慣性坐標系和入軌點軌道坐標系下的位置約束關系滿足

(29)

當入軌點軌道坐標系下的Y和Z方向的位置約束滿足，即E1=E2=0時，ΔXchi、ΔYchi和ΔZchi之間存在固定的比例關系，即無論(Xf(1)·Re-Xocff)取何值，均應有

(30)

此時地心慣性坐標系下等效的兩個終端約束為

(31)

(32)

在具體利用式(29)計算ΔXchi、ΔYchi和ΔZchi時，可直接將[Xf(1)·Re-Xocff]替換為某一常數C，同時，為進一步提高數值求解精度，對約束E1和E2可以進行相應的權重調整，即分配相應的權重系數k1和k2，從而有

(33)

類似地，對入軌點軌道坐標系下的3個速度約束E3、E4、E5也可以進行相應的權重調整，相應的權重因子取為k3、k4、k5。

總結來說，轉換后的5個等效終端約束為

(34)

在獲得等效的終端約束后，每一次制導循環中，通過對式(34)和式(17)、式(21)直接迭代求解來獲得協態變量初值，進一步由式(22)求解相應的控制角。

3　仿真驗證

3.1標準條件仿真

在發射慣性系下，航天器初始的位置為[1 865 014.8，40 816.2，150 433.5]Tm，速度為[7 412.601，-2 160.522，-130.991]Tm/s。Dfinal點的位置為[-5 967 060.6，-9 089 564.4，-89 871.6]Tm，速度為[-3 054.343，7 167.845，286.663]Tm/s。

仿真計算步長選為10ms，終止迭代計算條件選為剩余分析時間小于5s時，3次關機點條件為剩余飛行時間小于0.1s時。終端約束的權重因子取為k1=10-4、k2=10-4、k3=10-3、k4=10-3、k5=10-4。關機點的反饋系數k取為8。

首先對初選的開關機點直接進行制導仿真，制導結束后的偏差數據如表1所示。

表1　位置偏差和速度偏差(初始開關機點)

由表1中的數據可知，由于沒有對X方向的位置進行約束，因此最終X方向的位置偏差比較大，利用本文中的方法對關機點進行調整，迭代修正計算的次數為2次，第一次修正后X方向的位置偏差為25.806 7m，第二次修正后X方向的位置偏差為5.738 9m，可見反饋系數的選取能夠使得X方向的位置偏差朝減小的方向進行。

對修正后的開關機點進行制導仿真，實際飛行時間為166.890 0s，在入軌點軌道坐標系下的制導結果如圖3～6所示。

圖3　剩余時間示意Fig.3　Residual time sketch map

圖4　位置偏差Fig.4　Position deviation in the XYZ-direction

圖5　速度偏差Fig.5　Velocity deviation in the XYZ-direction

圖6　姿態角變化Fig.6　Attitude angleφandψ

由圖3～6的仿真結果可以看出，在整個制導過程中，剩余時間逐漸減小到0，3個方向的位置和速度誤差最終也趨于0，制導結束后的偏差數據如表2所示。

表2　位置偏差和速度偏差(修正開關機點)

從表1和表2可以看出，X方向的終端位置約束精度提高了159.535 5m,此后航天器開始自由滑行，滑行仿真在地心慣性坐標系下解算，滑行時間為745.1s，將最終到達Dfinal點時的仿真計算結果與所要求的相比結果如表3所示。

表3　到達Dfinal點時的偏差

由表3中的數據可知，航天器到達最終點時的各項偏差均在容許范圍內，最終的經度和緯度偏差均小于0.5°，高度偏差小于200m，速度傾角偏差小于0.005°，速度方向角偏差小于0.3°、速度大小偏差小于1.5m/s，從而表明了本文所提方法的有效性。

3.2蒙特卡羅打靶仿真

下面對航天器初始的位置和速度存在偏差的情況進行蒙特卡羅打靶仿真。假定位置和速度偏差相互獨立且服從正態分布，均值取為第3.1節中的數據，具體取值如下：

在發射慣性系下，航天器初始的位置為[1 865 014.8，40 816.2，150 433.5]Tm，速度為[7 412.601，-2 160.522，-130.991]Tm/s。

X方向位置偏差：均值μpx=1 865 014.8m，標準差σpx=1km。

Y方向位置偏差：均值μpy=40 816.2m，標準差σpy=100m。

Z方向位置偏差：均值μpz=150 433.5m，標準差σpz=300m。

X方向速度偏差：均值μvx=7 412.601m/s，標準差σvx=1m/s。

Y方向速度偏差：均值μvy=-2 160.522m/s，標準差σvy=0.5m/s。

Z方向速度偏差：均值μvz=-130.991m/s，標準差σvz=0.3m/s。

針對上述偏差數據，其余仿真設置與第3.1節相同，進行蒙特卡洛仿真，樣本點取為2000，某一次的仿真統計結果如圖7～9所示。

圖7　開關機點優化迭代次數Fig.7　Iteration times of switch on and off instants optimization

圖8　初始關機點X位置偏差大小Fig.8　Absolute position deviation of initial switch off instants in the X-direction

圖9　優化關機點X位置偏差大小Fig.9　Absolute position deviation of optimized switch off instants in the X-direction

進一步的統計結果如表4所示。

表4　制導結束時的偏差蒙特卡羅仿真統計結果

由圖7～9和表4中的統計結果可知，大部分的樣本點，經過關機點優化后能夠保證X方向位置偏差大小在40m以內，且迭代求解次數不超過5次。

最終到達Dfinal點時的仿真統計結果表5所示。

由表5中的統計結果可知，航天器到達最終點時的仿真，大約70%的樣本點能達到容許的偏差范圍，從而表明了所提方法對航天器初始位置和速度偏差具有一定的適應性。

表5　到達Dfinal點導的偏差蒙特卡羅仿真統計結果

4　結束語

本文針對傳統迭代制導方法的不足，首先對關機點進行優化，以減小未被滿足的終端位置約束的影響，進一步，直接從最優控制原理出發，推導控制角需要滿足的方程并進行求解，擺脫了傳統迭代制導小角度修正的假設，同時對終端約束進行等效轉換和權重分配，從而提高了數值求解的精度和制導方法的適應性，航天器最終能以指定精度到達下一次任務點，仿真結果表明了所提方法的有效性。

References)

[1]王召輝，金磊，賈英宏，等. 上面級軌道轉移段姿態控制技術[J]. 中國空間科學技術，2014，10(5): 41-48,55.

WANG Z H，JIN L，JIA Y H，et al. Attitude control technology of upper stage during orbital transfer[J]. Chinese Space Science and Technology，2014，10(5): 41-48,55(in Chinese).

[2]周洋，閆野，黃煦，等. 航天器軌道機動的閉環控制精度分析[J]. 宇航學報，2014，35(9): 1015-1021.

ZHOU Y，YAN Y，HUANG X，et al. Accuracy analysis of closed-Loop control for spacecraft orbit naneuver[J]. Journal of Astronautics，2014，35(9): 1015-1021(in Chinese).

[3]王常虹，曲耀斌，陸智俊，等. 航天器有限推力軌道轉移的軌跡優化方法[J]. 西南交通大學學報，2013，48(2): 390-394.

WANG C H，QU Y B，LU Z J，et al. Trajectory optimization method for spacecraft orbit transfer with finite thrust[J]. Journal of Southwest Jiaotong University，2013，48(2): 390-394(in Chinese).

[4]廖飛，季海波，解永春. 追蹤器本體坐標系下航天器姿軌一體化控制律設計[J]. 控制與決策，2015，30(9): 1679-1684.

LIAO F，JI H B，XIE Y C. Integrated orbit and attitude control for spacecraft in body fixed coordinate of chaser[J]. Control and Decision，2015，30(9):1679-1684(in Chinese).

[5]LU P，FORBES S，BALDWIN M. A Versatile powered guidance algorithm[C]∥2012 AIAA Guidance，Navigation，and Control Conference. Minnesota: Minneapolis，2012.

[6]HAEUSSERMANN W. Saturn launch vehicle′s navigation guidance，and control system[J]. Automatica，1971，7(5): 537-556.

[7]韓祝齋. 用于大型運載火箭的迭代制導方法[J]. 宇航學報，1983，4(1): 10-16.

HAN Z Z. An iterative guidance method for the large launch vehicle[J]. Journal of Astronautics，1983，4(1): 10-16(in Chinese).

[8]陳新民，余夢倫. 迭代制導在運載火箭上的應用研究[J]. 宇航學報，2003，24(5): 484-489.

CHEN X M，YU M L. Study of iterative guidance application to launch vehicles[J]. Journal of Astronautics，2003，24(5)： 484-489(in Chinese).

[9]茹家欣. 液體運載火箭的一種迭代制導方法[J]. 中國科學: E 輯，2009, 39(4): 696-706.

RU J X. An iterative guidance method of liquid launch vehicle[J]. Science in China Series E-Technological Sciences，2009 ,39(4): 696-706(in Chinese).

[10]LU P，GRIFFIN B J，DUKEMAN G A，et al. Rapid optimal multiburn ascent planning and guidance[J]. Journal of Guidance，Control，and Dynamics，2008，31(6): 1656-1664.

[11]LU P，SUN H，TSAI B. Closed-loop endoatmospheric ascent guidance[J]. Journal of Guidance，Control，and Dynamics，2003，26(2): 283-294.

[12]LU P，PAN B. Highly constrained optimal launch ascent guidance[J]. Journal of guidance，Control，and Dynamics，2010，33(2): 404-414.

[13]吳楠，程文科，王華. 運載火箭迭代制導方法的改進研究[J]. 動力學與控制學報，2009，7(3): 279-282.

WU N，CHENG W K，WANG H. An improved iterative guidance method for launch vehicle[J]. Journal of Dynamics and Control，2009，7(3): 279-282(in Chinese).

[14]王煬，唐碩. 最優真空飛行軌跡設計和制導方法研究[J]. 計算機仿真，2011，28(10): 78-82.

WANG Y，TANG S. Optimal vacuum flight trajectory design and guidance method[J]. Computer Simulation，2011，28(10):78-82(in Chinese).

[15]傅瑜，陳功，盧寶剛，等. 基于最優解析解的運載火箭大氣層外自適應迭代制導方法[J]. 航空學報，2011，32(9): 1696-1704.

FU Y，CHEN G，LU B G，et al. A vacuum adaptive iterative guidance method of launch vehicle based on optimal analytical solution[J]. Acta Aeronautica et Astronautica Sinica，2011，32(9): 1696-1704(in Chinese).

(編輯：車曉玲)

An optimal guidance method based on multiple terminal constraints

LI Chaobing*，WANG Jinlin，LI Hai

Beijing Aerospace Automatic Control Institute，Beijing 100854，China

Since the thrust amplitude of the main engine on spacecraft is unadjustable，a multiple terminal constrain based optimal guidance method was proposed，considering the deficiency of small correction angle assumption of conventional iterative guidance under five terminal constrains. An optimal control model was established in the orbit injection orbital coordinate frame. Transversality equations ware directly solved by iteration to obtain the guidance angle command.The switch on and off instants were then optimized to reduce the influence of unsatisfied terminal constraints. Besides，the five equivalent terminal constraints in the geocentric inertial coordinate frame were derived and appropriate weight factors were adjusted to improve the precision of numerical solving of guidance equations. Simulations on standard conditions demonstrate that the proposed guidance method has increased the precision of unsatisfied terminal constraints by 159.535 5 m，compared with the conventional iterative guidance. The Monte Carlo simulations show the adaptability of the proposed method to the initial bias of the spacecraft position and velocity.

spacecraft orbit transfer；iterative guidance；optimal control；coordinate transformation；switch on and off instants optimization；weight adjustment

10.16708/j.cnki.1000-758X.2016.0052

2016-05-03；

2016-07-07；錄用日期：2016-08-22；

時間：2016-09-2113:41:23

http:∥www.cnki.net/kcms/detail/11.1859.V.20160921.1341.005.html

李超兵(1981-)，男，碩士，高級工程師，lcbpku@163.com，研究方向導航制導與控制

TP249，V448.22

A

http:∥zgkj.cast.cn

引用格式：李超兵，王晉麟，李海. 一種基于多終端約束的最優制導方法[J].中國空間科學技術，2016，36(5)：9-17.

LICB，WANGJL，LIH.Anoptimalguidancemethodbasedonmultipleterminalconstraints[J].ChineseSpaceScienceandTechnology，2016，36(5)：9-17(inChinese).