脈沖推力最優軌跡的Hamilton邊值問題

沈紅新

(西安衛星測控中心宇航動力學國家重點實驗室，西安 710043)

?

脈沖推力最優軌跡的Hamilton邊值問題

沈紅新

(西安衛星測控中心宇航動力學國家重點實驗室，西安 710043)

針對大推力航天器的Hamilton邊值問題(HBVP),提出一組基于變分法的通用方程，其中內點和其它端點(包括始、末端點)可以滿足統一的方程形式，由此反映了更本質的邊值條件解析結構。具體問題的最優性必要條件均可以從本文給出的通用方程中較方便地推出，避免了以往構造邊值問題復雜繁瑣的困難。仿真結果表明，本文方法可以保證有效、快速地獲得大推力航天器的最優飛行路徑。

邊值問題; 最優性；內點；脈沖軌道

0 引 言

航天任務的經濟、技術可行性和航天器飛行路徑密切相關，最小化燃料和時間資源的消耗是航天器設計者不懈的追求，因此有效的優化控制方法對任務設計和方案制定至關重要。很多航天任務特別是載人任務都采用大推力的化學推進，為此，本文研究大推力航天器的最優飛行路徑問題。目前，在非線性最優控制領域主要有兩大類數值方法：直接法和間接法[1-3]。直接法一般采用參數化方法將連續空間的非線性最優控制問題轉換成離散空間的非線性規劃問題，雖然它具有容易實現和收斂半徑較大的優點，但是直接法不能保證所獲得解的最優性。相對地，間接法主要是構造和求解由最優控制的一階必要條件得到的Hamilton邊值問題(Hamilton boundary value problem,HBVP)[3-4]，其優點是：1)構造邊值問題而不是參數優化問題，理論上能夠準確快速獲得最優解；2)由于采用最優控制理論(這里主要指變分法和極大值原理)，間接法能夠提供優化問題必要的理論信息。但間接法也有其缺點：獲得最優性條件的過程復雜繁瑣，而且收斂半徑較小[3]。

最優飛行路徑確定的核心問題是構造和求解HBVP。近年來有很多研究集中于對HBVP的求解并產生了一批研究成果，如辛方法[5]、生成函數方法[6]、同倫變換[7-9]和微分變換[10]等，但這些研究都較少關注HBVP的構造方法。針對HBVP確立最優性條件復雜繁瑣的缺點，本文的目的是以比較簡潔的方程形式給出大推力航天器最優路徑規劃的通用模型，為該類問題最優性條件的建立提供統一的理論框架，這將有助于設計者快速、方便地構造進而求解HBVP。在任務初步設計階段，大推力一般被假設為瞬時脈沖作用，因此大推力航天器最優路徑確定是典型的含有內點約束的最優控制問題。Bryson和Ho[4]以及Lawden[11]研究了狀態不連續內點約束最優控制問題，但他們的研究僅給出了內點處的最優性條件。本文將對內點和其它一般端點(例如初始時刻、終端時刻等)給出統一的邊值問題模型，這樣做的好處是得到的邊界條件更具有一般性，可以方便確定起始時刻t0+和終端時刻tf-，而這兩個時刻在軌跡優化問題中也可能是未知量 (例如自由終端時刻的問題)；相對地，Bryson和Ho[4]以及Lawden[11]的理論只適用于固定時間軌道轉移或交會問題。

1 問題提出

由于化學推進的推力很大，發動機工作時間相對飛行時間來說很短甚至可以忽略，所以可假設化學推進的推力為瞬時脈沖，脈沖作用前后航天器位置連續而速度發生跳變，脈沖作用點之間由無動力滑行軌跡聯接。無動力滑行軌跡滿足如下動力學系統微分方程：

(1)

式中：x為狀態變量，t為自變量(本文中自變量為時間)。需要指出的是，這里脈沖并不作為控制變量，而只是作為狀態間斷來考慮，所以在脈沖軌跡最優控制中不存在控制變量。

在軌跡初始點和終點之間如果存在狀態或控制不連續，或者存在約束，這樣的位置被稱為內點，在內點處需要給出對應的最優性必要條件。比較方便的做法是將軌跡根據內點劃分為多段，例如第j段從t(j-1)+開始到tj-結束，對應的狀態分別是x(j-1)+和xj-，其中j-和j+分別代表在點j之前和之后。脈沖軌跡由于存在速度狀態的突變，所以是典型的含有內點約束問題。在軌跡始末兩個端點和內點施加非線性約束條件，這些約束用一個函數ψ描述，其表達式寫為

ψ(x(j-1)+,xj-,t(j-1)+,tj-)=0,j=1,…,f

(2)

從容許函數類中求一函數x(t)，使泛函

J=φ(x(j-1)+,xj-,t(j-1)+,tj-),j=1,…,f

(3)

取極大(小)的問題常稱為Mayer問題，它是古典變分的三個基本問題之一。此外還有Bolza問題和Lagrange問題。引進某些輔助變量可使三類問題互相轉化。因此，研究三種基本問題中的任何一種都具有普遍意義，本文采用Mayer形式來定義優化問題。此外，通過改變目標函數的數學形式，求泛函極小的問題可以轉化為求泛函極大的問題，而本文只考慮求極大值的問題。

對不同優化指標及約束條件的問題，HBVP確立最優性條件比較復雜繁瑣[3]，本文的目的是為脈沖轉移這類問題最優性條件的建立提供簡潔的統一的理論框架，這將有助于設計者快速、方便地構造進而求解HBVP。

2 Hamilton邊值問題構造的一般形式

變分法一般從推演泛函極值的必要條件開始。 因為式(3)表示的性能指標與約束條件或微分方程無關， 故引入和邊界條件相聯系的Lagrange乘子μ， 以及和微分方程相聯系的協變量λ， 從而構造新的泛函J*：

(4)

(5)

從物理意義的角度說，只用時間和狀態對這個泛函進行微小擾動，如果擾動的結果是這個泛函不動，說明是泛函達到了極值，也就說明達到了最佳的時間和狀態。當然如果泛函中包含控制，必然會有控制的變分，但本文中脈沖作用被視為狀態的間斷而非控制量。

分部積分式(5)最后一項得

(6)

定義Hamilton函數為

H=λTf

(7)

并考慮到狀態量的變分和自變量有關，有

(8)

這里認為自變量t的微分等于變分。

將式(6)～(8)代入式(5)得

(9)

最優性必要條件應該由極值的基本必要條件確定，即要求泛函J(或J的等價泛函J*)的一階變分為0，即

δJ*=0

(10)

并且不依賴于dtj-、dt(j-1)+、dxj-、dx(j-1)+和δx的選擇。通過使δx的系數為0，獲得Euler-Lagrange方程

(11)

Euler-Lagrange方程是描述協變量演化的微分方程，和狀態方程相對應，也稱為協變量方程。

其它系數和每段端點有關，令dtj-、dt(j-1)+、dxj-、dx(j-1)+的系數為0得到一般性的最優邊界條件為

(12)

(13)

(14)

(15)

式(12)、式(13)和確定最優狀態有關，式(14)、式(15)和確定最優時間有關。式(14)和式(15)在t(j-1)+和tj-時刻分別得到了Hamilton函數的兩個邊界條件。

與Bryson給出的關于內點時刻tj的一個條件相比，在討論起始時刻t0+和終端時刻tf-時，式(14)和式(15)更為明確，因為這樣避免了涉及t0-和tf+這些實際不存在的時刻點。內點(即脈沖作用點)在上述方程中只是作為一個普通的端點處理，內點和其它端點(始、末端點)可以滿足統一的方程形式。

內點的一階必要條件和Lawden條件相同，但Lawden必要條件沒有給出端點(初始時刻和終端時刻)一階必要條件，而在間接法中，端點和內點的必要條件可以從同一組方程中導出。相對Bryson和Ho[4]以及Lawden[11]給出的內點方程而言，本文的Hamilton邊值問題方程更通用，由此也反映了更本質的邊值條件解析結構。具體問題的最優性必要條件均可以從本文給出的通用方程中較方便地推出，避免了以往構造邊值問題復雜繁瑣的困難。

需要說明的是，本文構造的含內點約束的多點邊值問題(Multi-pointboundaryvalueproblem，MPBVP)總能得到和未知量相同數目的邊界條件，具體原因分析如下：

1)式(12)～(15)中μ包含的每個Lagrange乘子都為常數，由于μ和約束ψ對應，所以它的個數等于ψ所包含約束的個數。問題求解過程中，為了減少未知量的個數，一般會通過一些代數手段消掉Lagrange乘子。

2)式(12)、(13)能夠為初始狀態x0+和初始協變量λ0+提供數目相等的邊界條件(j分別等于f和0)。

3)式(14)、(15)能夠為t(j-1)+和tj-提供數目相等的邊界條件。

4)式(12)、(13)能夠為內點未知量提供數目相等的邊界條件。

3 仿真校驗

本節采用三個典型算例來測試上文給出的HBVP通用模型，分別是霍曼轉移、異面變軌和多脈沖交會。

由于球坐標形式比直角坐標形式物理概念更為直觀，而且下文也會展示采用球坐標形式時，其中一個協變量始終為常值，因此本文算例采用球坐標運動方程。只考慮地球引力，航天器的運動方程和對應的協變量微分方程為[12]

(16)

(17)

(18)

(19)

(20)

(21)

λvθ(-vrvθ+vθvφtanφ)+λφv+

(22)

(23)

(24)

(25)

(26)

(27)

3.1 霍曼轉移

第一個例子是共面圓軌道之間的軌道轉移問題。考慮初始軌道為半徑r0=8000 km的環火星圓軌道，目標軌道為半徑rf=15000 km的圓軌道。這個問題同Abdelkhalik和Mortari[13]采用遺傳算法結合最速下降法計算脈沖軌道轉移的第一個算例相同，該算例源自文獻[14]。Brown[14]解析證明了此問題最優解的總速度增量為0.609km/s，轉移時間為5.2h。

目標函數為總的速度增量最小，最大化性能指標寫為

(28)

式中：下標“0”、“2”和“f”分別表示第一次脈沖之前、第二次脈沖時間之前和之后的狀態。初始和終端條件分別為

(29)

r2=rf

(30)

根據式(13)和式(28)得邊值條件

(31)

(32)

根據式(12)和式(28)得邊值條件

(33)

(34)

由于終端時刻自由，根據式(14)得Hamilton函數邊值條件

H2=0

(35)

由式(29)和式(13)得，λr1=μr1，其中μr1是對應等式約束r1=r0的Lagrange乘子，所以λr1是未知參數。θ0的取值是任意的，不妨設θ0=0；由于終端條件不限制θ，有λθ2=0，結合式(23)可知，λθ是一個等于0的常數。綜上，式(29)～(35)給出了8個邊界條件。優化模型的設計變量為8個，分別為Δt、r1、θ1、vr1、vθ1、λr1、λvr1和λvθ1，其中下標“1”表示第一次脈沖之后的狀態，未知量和設計變量的數目相等，上述優化問題被轉化為兩點邊值問題。

本文中邊值問題求解方法采用解析同倫方法，采用這種方法可以從具有解析解的初始構造問題出發，通過調整參數逐步迭代過渡到原始優化問題，因此這種方法的優勢是可以避免協變量初值的猜測的困難，從而比較容易求解邊值問題。由于邊值問題求解不是本文的研究重點，方法的原理參見文獻[12]，此處不再贅述。最終得到協變量初值為λr1=0.268，λvr1=0，λvθ1=1，計算所得結果如表1所示，和理論上的霍曼轉移解析解相同。Abdelkhalik和Mortari[13]報告了一個接近霍曼轉移的次優解，其中所需速度增量和時間分別為0.6093km/s和5.157h，和最優解的偏差分別為0.033%和0.83%。需要指出，Abdelkhalik和Mortari[13]所計算的霍曼轉移最優飛行時間5.08h有誤，應為5.2h，見表1。

圖1描述了終端時間自由時的最優解主矢量變化歷史，顯然滿足主矢量必要條件。可見，采用間接法能夠得到滿足主矢量條件的最優解。

表1 霍曼轉移問題求解結果

圖1 霍曼轉移問題最優解的主矢量變化歷史Fig.1 Primer history of the optimal LEO to HEO transfer

3.2 異面變軌

第二個例子是異面圓軌道之間的軌道轉移問題。 考慮初始軌道為半徑為6671.53km的地球停泊圓軌道，軌道傾角28.5°，目標軌道為半徑26558.56km的圓軌道，軌道傾角為0°。這個問題同Abdelkhalik和Mortari[13]采用遺傳算法結合最速下降法計算脈沖軌道轉移的第二個算例相同，該算例源自文獻[15]。Vallado[15]提供了兩脈沖最優異面變軌的解析解法。

目標函數為總的速度增量最小，最大化性能指標寫為

(36)

式中：下標“0”、“2”和“f”分別表示第一次脈沖之前、第二次脈沖時間之前和之后的狀態。

考慮這個問題的一種限制性情況，即限制轉移航天器初始時刻在兩個軌道的交線上，則φ0=0、vr0=0、vθ0=vc1cosi、vφ0=vc1sini，其中vc1和i分別表示初始軌道的速度大小和軌道傾角，不妨設θ0=0。初始和終端條件分別為

(37)

(38)

根據式(13)和式(36)得邊值條件

(39)

(40)

(41)

根據式(12)和式(36)得邊值條件

(42)

(43)

(44)

和前面的共面轉移問題類似，有

H2=0

(45)

且λθ恒等于0，

(46)

式中：μr1和μφ1是引入的Lagrange乘子，是未知參數。

綜上，式(37)～(45)給出了12個邊值條件。優化模型的設計變量分別為Δt、r1、θ1、φ1、vr1、vθ1、vφ1、λr1、λφ1、λvr1、λvθ1和λvφ1，其中下標“1”表示第一次脈沖之后的狀態。可見，設計變量的個數和邊界條件的個數相等，上述優化問題被轉化為兩點邊值問題。

邊值問題求解方法同樣采用解析同倫方法，得到的協變量初值為λr1=0.757，λφ1=6.16×10-8，λvr1=-6.84×10-8，λvθ1=0.975，λvφ1=0.224，計算結果如表2所示。本文方法能夠得到和理論最優解相同的結果，并且優于采用遺傳算法結合最速下降法獲得的結果(總的速度增量為4.0610km/s[13])和采用偽譜法獲得的結果(總的速度增量為4.0594km/s[16])，見表2。圖2描述了終端時間自由時的最優解主矢量變化歷史，顯然滿足主矢量必要條件。

表2 異面變軌問題求解結果

圖2 異面變軌問題最優解的主矢量變化歷史Fig.2 Primer history of the optimal noncoplanar transfer

3.3 多脈沖交會

同圓交會問題是一個比較經典的多圈脈沖交會問題。目標和追蹤航天器是在同一個圓軌道上，追蹤航天器滯后目標航天器一定的相位角，要求在一定時間內追蹤航天器與目標航天器在該圓軌道上交會(相對狀態為0)。圖3是目標航天器在追蹤航天器前方180°的示意圖。關于這個交會問題的研究已經很多，包括Prussing和Chiu[17]、Colasurdo等[18]、Prussing[19]以及Luo等[20-21]，當交會時間是2.3個軌道周期時，上述研究分別報告了最優四脈沖交會解，而且該交會問題需要一個非常大的Δv(大于1200 m/s)。為了便于對比，本文也選擇一個400 km圓軌道工況進行測試。

多脈沖交會問題的HBVP最優性條件推導方法和前兩個問題類似，因此推導過程這里不再重復，而是直接給出設計結果，特別是主矢量變化結果，從中可以觀察解的最優性必要條件的滿足情況。

第一個四脈沖最優解主矢量變化見圖4(圖中黑點表示脈沖作用點)，顯然滿足Lawden條件。需要指出，這個四脈沖最優解雖然滿足主矢量條件，但并不是可行解，因為航天器飛行軌跡最小距離小于地球半徑(見圖5)。

第二個四脈沖最優解主矢量變化見圖6(圖中黑點表示脈沖點)，顯然滿足Lawden條件。這組四脈沖解不但滿足主矢量條件，而且也是可行解，圖7給出了航天器飛行軌跡在赤道面內的投影，可見最小距離大于地球半徑。

圖3 追蹤和目標航天器在同一圓軌道上交會問題Fig.3 Illustration of the same-circle rendezvous problem

圖4 同圓交會問題第一個解主矢量(算列1)Fig.4 Primer history of same-circle rendezvous (case 1)

圖5 同圓交會問題第一個四脈沖解最優軌跡(算列1)Fig.5 Optimal solution of same-circle rendezvous (case 1)

圖6 同圓交會問題第二個解主矢量(算列2)Fig.6 Primer history of same-circle rendezvous (case 2)

圖7 同圓交會問題第二個四脈沖解最優軌跡(算列2)Fig.7 Optimal solution of same-circle rendezvous (case 2)

Prussing和Chiu[17]最早對同圓交會問題進行了研究，但只給出了第二個四脈沖最優解。Colasurdo等[18]采用間接法并結合動力學輔助分析首次完整給出了上述兩個四脈沖最優解。Luo等[20]和Luo等[21]利用進化算法也找到了上述兩個四脈沖最優解，然而進化算法通常需要的計算較大。相比之下，本文提出的方法較為簡單通用，而且能從理論上保證解的最優性。

需要說明的是，初始和終端時刻不能作為內點處理，而且終端時刻可能不固定(見前兩個算例)，因此主矢量理論不適用于這兩個端點，本文構造的模型可以將端點和內點統一起來處理。

以上三個算例都比較典型，代表了不同類型的軌道轉移問題：前兩個算例有解析解，一個是平面轉移，另一個是異面轉移，第三個算例是典型的多脈沖交會問題，不存在解析解。用不同類型的問題能夠更好地驗證文中模型的適用性。其他脈沖轉移問題的Hamilton邊值問題的構建與此類似。Lawden主矢量理論是在二體模型中推導的，本文提出的Hamilton邊值問題通用模型并不依賴于二體動力學模型假設，對更加復雜的動力學問題也有應用前景。

表3 400 km同圓交會軌道問題的二組解

4 結 論

本文研究了大推力航天器最優飛行路徑的Hamilton邊值問題，給出了這類問題所滿足的邊值條件一般形式。根據邊值條件的一般形式，可以很方便地得到具體問題的邊值條件，從而簡化了原問題的復雜性。在討論最優飛行路徑Hamilton邊值問題的過程中，發現內點和所有其它端點均可以滿足同樣一套方程，據此，拓展和完善了Bryson和Ho求解該問題的模型。由于采用了脈沖假設，這個假設只適用于大推力情形，對于小推力問題，是否存在比較通用的邊值條件形式還值得進一步研究。

[1] 雍恩米, 陳磊, 唐國金. 飛行器軌跡優化數值方法綜述[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.]

[2] 李俊峰, 蔣方華. 連續小推力航天器的深空探測軌道優化方法綜述[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.]

[3] Betts J T. Survey of numerical methods for trajectory optimization [J]. Journal of Guidance, Control, and Dynamics, 1998, 21(2): 193-207.

[4] Bryson A E, Ho Y C. Applied optimal control[M]. Washington, DC: Hemisphere, 1975.

[5] Peng H J, Gao Q, Wu Z, et al. Symplectic approaches for solving two-point boundary-value problems [J]. Journal of Guidance Control and Dynamics, 2012, 35 (2): 653-658.

[6] 陳琪鋒,張躍東,吳文昭,等. 基于近似生成函數迭代的分布式衛星構形最優控制研究[J]. 宇航學報, 2009, 30(3): 988-993. [Chen Qi-feng, Zhang Yue-dong, Wu Wen-zhao, et al. Research on optimal control of distributed satellite system formation based on iteration of generating function approximation [J]. Journal of Astronautics, 2009, 30(3): 988-993.]

[7] Jiang F, Baoyin H, Li J. Practical techniques for low-thrust trajectory optimization with homotopic approach [J]. Journal of Guidance, Control, and Dynamics, 2012, 35(1): 245-258.

[8] Shen H X, Casalino L, Li H Y. Adjoints estimation methods for impulsive Moon-to-Earth trajectories in the restricted three-body problem [J]. Optimal Control Applications and Methods, 2014, 36: 463-474.

[9] Pan B, Lu P, Pan X, et al. Double-homotopy method for solving optimal control problems [J]. Journal of Dynamic Systems, Measurement and Control, 2016, 39(8): 1706-1720.

[10] Hwang I, Li J H, Du D. Differential transformation and its application to nonlinear optimal control [J]. Journal of Dynamic Systems, Measurement and Control, 2009, 131(5): 1-11.

[11] Lawden D F. Optimal trajectories for space navigation [M]. London: Butterworths, 1963.

[12] 沈紅新. 基于解析同倫的月地應急返回軌跡優化方法[D]. 長沙：國防科技大學, 2014. [Shen Hong-xin. Optimization method for the Moon-Earth abort return trajectories based on analytic homotopic technique [D]. Changsha: National University of Defense Technology, 2014.]

[13] Abdelkhalik O, Mortari D. N-impulse orbit transfer using genetic algorithms [J]. Journal of Spacecraft and Rockets, 2007, 44(2): 456-460.

[14] Brown C D. Spacecraft mission design[M]. AIAA, 1998.

[15] Vallado D A. Fundamentals of astrodynamics and applications[M]. Torrance, USA: Microcosm Press, 2001.

[16] Shen H X. A novel algorithm for optimal design of impulsive orbit transfer [C]. The 3rd International Symposium on Systems and Control in Aeronautics and Astronautics, Harbin, China, June 8-10, 2011.

[17] Prussing J E, Chiu J H. Optimal multiple-impulse time-fixed rendezvous between circular orbits [J]. Journal of Guidance Control and Dynamics, 1986, 9(1):17-22.

[18] Colasurdo G, Pastrone D. Indirect optimization method for impulsive transfer[R]. AIAA Paper 1994-3762, 1994.

[19] Prussing J E. A class of optimal two-impulse rendezvous using multiple-revolution Lambert solutions [J]. The Journal of the Astronautical Sciences, 2000, 48(2):131-148.

[20] Luo Y Z, Zhang J, Li H Y, et al. Interactive optimization approach for optimal impulsive rendezvous using primer vector and evolutionary algorithms[J]. Acta Astronautica, 2010, 67(3): 396-405.

[21] Luo Y Z, Tang G J, Li Y J, et al. Optimization of multiple-impulse, multiple-revolution, rendezvous phasing maneuvers [J]. Journal of Guidance, Control, and Dynamics, 2007, 30(4): 946-952.

[22] Shen H X, Casalino L. Indirect optimization of three-dimensional multiple-impulse Moon-to-Earth transfers [J]. The Journal of Astronautical Sciences, 2014, 61: 255-274.

通信地址：西安505信箱28分箱(710043)

電話：(029)84762520

E-mail:shxnudt@163.com

Hamilton Boundary Value Problem for Optimal Impulsive Trajectory

SHEN Hong-xin

(State Key Laboratory of Astronautic Dynamics, Xi’an Satellite Control Center, Xi’an 710043, China)

For a spacecraft with high thrust, a general equations based on variational calculus for formulating the Hamilton boundary value problem (HBVP) is proposed, where the optimality conditions on interior points and other general boundary points (inclucling initial and final points) are in the same manners, so that our equations tend to reveal the more intrinsic analytic structure of HBVP. The detailed optimality conditions can be determined conveniently according to the general equations; therefore, the tedious task of formulating HBVP is precluded. Simulations show that the method proposed in this paper could ensure obtaining the optimal flight path of the high-thrust spacecraft effectively and efficiently.

Boundary value problem; Optimality; Interior point; Impulsive trajectory

2016-10-18;

2017-05-15

V412.4

A

1000-1328(2017)07-0686-08

10.3873/j.issn.1000-1328.2017.07.000

沈紅新(1986-)，男，博士，工程師，主要從事航天器軌道姿態最優控制研究。