袁 晗，王小軍，牟 宇，張宏劍，張普卓

(1. 北京宇航系統工程研究所，北京 100076；2. 中國運載火箭技術研究院，北京 100076)

0 引 言

航天運載器的重復使用具有降低發(fā)射成本、解決落區(qū)安全的重要意義，2015年以來美國的私營航天公司SpaceX和Blue Origin已經通過垂直起降的方式成功實現火箭一子級回收和重復使用，驗證了垂直起降重復使用火箭(簡稱垂直起降火箭)在目前的技術水平下是一種技術上可行、經濟上合理的運載器重復使用方案。設計垂直起降火箭返回的制導控制方案時，通常將飛行過程劃分為多個階段分別進行設計。以不返回發(fā)射場的任務為例，其返回飛行剖面如圖1所示，各階段分別為：調姿滑行段AB、動力減速段BC、氣動減速段CD和動力著陸段DE等。其中，動力著陸段具有初始狀態(tài)偏差大、著陸精度要求高，且燃料消耗受制導方法影響大等特點，近年來諸多學者對動力著陸制導開展了大量的研究，通常基于在線軌跡規(guī)劃，通過間接法或直接法等軌跡規(guī)劃方法在線生成當前狀態(tài)對應的最優(yōu)軌跡，從而獲得制導指令輸出。為降低重復使用帶來的火箭運載能力損失，制導方法中通常以燃料最優(yōu)作為優(yōu)化目標。燃料消耗不僅受制導方法影響，還與動力著陸段的初始狀態(tài)有較大關系，不合理的啟動點將增加燃料消耗，甚至導致啟動點的狀態(tài)在動力著陸可行域之外。目前，部分研究通過單一狀態(tài)量觸發(fā)方法在線確定動力著陸段的啟動點，如時序觸發(fā)、高度觸發(fā)等。然而，初始狀態(tài)偏差較大時，上述方法難以保證啟動點為燃料最優(yōu)，且存在啟動點接近可行域邊界，甚至在可行域之外的風險。

圖1 典型垂直起降火箭飛行剖面Fig.1 Typical flight profile of vertical-takeoff-vertical-landing rockets

為解決單一狀態(tài)量觸發(fā)方法的不足，眾多學者對動力著陸段自適應啟動方法進行了研究。針對火星著陸問題，Lu等提出軟著陸軌跡預測自適應啟動方法，該方法通過在線預測軟著陸軌跡的航程判斷是否啟動發(fā)動機并啟動動力著陸段(簡稱啟動發(fā)動機)，軟著陸軌跡為僅施加終端高度和速度約束、松弛終端水平位置約束時的燃料最優(yōu)軌跡。但燃料最優(yōu)啟動點的動力著陸軌跡并非軟著陸軌跡，且二者航程存在顯著差異，故該方法難以得到燃料最優(yōu)啟動點。Jiang等提出強化學習決策方法，離線訓練決策是否啟動發(fā)動機的神經網絡策略函數，但神經網絡策略函數的可解釋性較差，導致目前工程上難以應用。針對火箭垂直回收問題，Wang等和邵楠等提出多種數值優(yōu)化方法，通過對氣動減速段和動力著陸段聯合軌跡規(guī)劃尋找燃料最優(yōu)啟動點，但由于在線計算效率難以滿足工程要求和算法收斂性難以保證，限制了該類方法的工程應用。Song等提出中間推力軌跡自適應啟動方法并通過GNC驗證飛行器“Peacock”進行了實驗驗證，該方法通過判斷火箭能否以恒定中間推力的軌跡實現著陸判斷是否啟動發(fā)動機，可使啟動點具有較高的魯棒性，但難以保證啟動點的燃料最優(yōu)性。

針對氣動減速段中在線判斷燃料最優(yōu)啟動點的問題，本文對火箭垂直返回中的燃料最優(yōu)啟動點進行了研究，采用一種解析的推力剖面和攻角剖面形式描述從燃料最優(yōu)啟動點出發(fā)的燃料最優(yōu)動力著陸軌跡(簡稱最優(yōu)啟動點的動力著陸軌跡)。結合上述剖面形式，提出了燃料最優(yōu)啟動點的判斷條件，進而設計了一種動力著陸段自適應啟動方法。仿真結果表明，該方法得到的啟動點與數值優(yōu)化方法得到的燃料最優(yōu)啟動點十分接近。

1 火箭著陸問題描述

針對研究的氣動減速段和動力著陸段的制導問題，本章將給出火箭返回的動力學模型，及氣動減速制導和動力減速制導的方案。

1.1 火箭返回的動力學模型

在以目標著陸點為原點的站心直角坐標系(簡稱站心坐標系)研究返回過程，該坐標系的向為縱向，向為高度，向為橫向。視地球為球形，忽略引力攝動。采用三自由度質點動力學模型，火箭的控制量為推力、攻角和側滑角，假定推力方向始終沿箭體軸線。動力學方程為

(1)

式中：為火箭的飛行速度；為速度傾角；為航跡偏航角；為火箭質量；為重力加速度；為火箭發(fā)動機比沖；,和為地球曲率、牽連加速度和Coriolis加速度在三個運動方向分別的影響，具體見文獻[17]；=(+)為重力加速度，其中為目標著陸點處的重力加速度，為目標著陸點處的地心距離；,和分別為氣動阻力、升力和橫向力，將其表示為

(2)

式中：為動壓；為阻力系數；′為升力系數導數，設為常數，升力系數()形式較為復雜難以解析表示，小攻角下可近似為()≈′，由于火箭為軸對稱體，橫向力系數()≈′；為大氣密度，采用標準大氣參數；為參考面積。

1.2 約束和優(yōu)化目標

以氣動減速段開始作為分析的初始狀態(tài)，設該時刻為=0，動力著陸段初始時刻為，終端時刻為。記=0時刻的狀態(tài)為[,,,,,,]。動力著陸段的初始狀態(tài)為氣動減速段的終端狀態(tài)。動力著陸段的終端位置、速度、傾角固定，終端質量和航跡偏航角無約束，動力著陸段的終端約束為

(3)

式中:≥0。以燃料最優(yōu)為目標，目標方程為

min=-()

(4)

氣動減速段無發(fā)動機推力，即=0。動力著陸段中發(fā)動機的最大推力、最小推力和推力調節(jié)速度均面臨約束，將其簡化為常值約束，分別為

≤≤

(5)

(6)

其中,>。將氣動載荷、舵面配平能力等假定為攻角和側滑角及其變化率的常值約束，為

(7)

(8)

終端時刻推力的大小無約束，由于終端時刻要求以豎直姿態(tài)觸地，故終端時刻的攻角、側滑角約束為

()=()=0

(9)

1.3 分段制導方案

火箭在氣動減速段中通過該段的制導算法生成攻角、側滑角指令，同時通過自適應啟動算法自主地確定是否啟動發(fā)動機。啟動發(fā)動機后，進入動力著陸段，且不再執(zhí)行自適應啟動算法。下面給出分析動力著陸啟動問題時采用的分段制導方案，其中氣動減速制導用于生成定量的返回軌跡，動力著陸制導用于生成動力著陸軌跡，進而給出啟動點的可行性和燃料消耗。

氣動減速制導采用文獻[17]的預測制導并進行改進。文獻[17]給定虛擬目標點在站心坐標系中坐標為[,,]，通過參數和確定預測軌跡的攻角剖面和側滑角剖面，利用對動力學方程式(1)積分預測=時的水平位置坐標|=,|=，進而以|==,|==為條件求解和。文獻[17]中設定攻角、側滑角為時間的線性函數，由于飛行時間受初始狀態(tài)和攻角、側滑角剖面影響較大，該方法限制了水平機動能力。因此不同于文獻[17]，本文設定攻角和側滑角指令為高度的函數。預測落點坐標時，式(1)中推力取0，攻角和側滑角與和的關系分別為

=min{max{(-)|-|,0},}

(10)

=min{max{(-)|-|,0},}

(11)

式中:,為預先給定的參數，即<時指令為0，<時指令為0，取=1 km,=5 km。

本文提出的自適應啟動方法不依賴于具體的動力著陸制導算法，采用數值軌跡優(yōu)化方法直接得到從啟動點出發(fā)的燃料最優(yōu)動力著陸軌跡來評價啟動點的可行性和燃料消耗，避免了不同制導算法對啟動點的可行性和燃料消耗評價結果的影響。軌跡優(yōu)化采用GPOPS算法，動力學模型為式(1)，終端約束為式(3)，優(yōu)化目標為式(4)。

2 燃料最優(yōu)啟動點的性質分析

在縱向平面(--平面)分析典型返回軌跡中使動力著陸段燃料消耗最小的啟動點(即燃料最優(yōu)啟動點，簡稱最優(yōu)啟動點)，及對應的動力著陸軌跡的性質。為便于得到控制量剖面的規(guī)律，文中忽略推力變化率約束(6)、攻角變化率約束(8)，終端狀態(tài)約束(3)中忽略傾角約束，僅考慮如下終端約束：()=0,()=0,()=,()=0。但在仿真校驗中考慮了上述忽略的約束，驗證了依據本章結論提出的算法仍有效。

采用如表1所示的火箭參數和初始狀態(tài)，通過如下方式生成多條氣動減速段軌跡：初始狀態(tài)相同，通過為氣動減速段的預測制導選取不同虛擬目標點生成不同軌跡，令虛擬目標點的高度均為=0。由于氣動減速制導采用無動力軌跡預測落點，無動力軌跡航程大于動力著陸航程，故虛擬目標點應在目標著陸點前方。數值仿真表明<786m時氣動減速段不存在動力著陸可行解，故取分別為100 m，150 m，…，500 m，共9條軌跡。采用GPOPS計算從氣動減速軌跡的不同時刻啟動對應的燃料消耗，采用啟動點的高度代表該啟動點，得燃料消耗與啟動高度曲線如圖2所示，圖中每條曲線下方為該曲線對應取值，由左至右依次對應為100 m，150 m，…，500 m的軌跡。

表1 火箭參數和初始狀態(tài)Table 1 Rocket parameters and initial states

圖2 不同虛擬目標點下燃料消耗與啟動點關系圖Fig.2 Relation between fuel consumption and starting points under different virtual target points

由圖2可見，最優(yōu)啟動點存在位于和不位于最小可行高度兩種情況：若位于最小可行高度，返回軌跡中隨啟動高度降低燃料消耗單調降低，如圖2中為100 m至200 m的3條軌跡；若不位于最小可行高度，返回軌跡中隨高度降低燃料消耗先降低后上升，如圖2中為250 m至500 m的6條軌跡。以分別為200 m，300 m和400 m的三條軌跡(圖2中三條加粗黑線)為例分析最優(yōu)啟動點的動力著陸軌跡。上述三個典型返回軌跡的最優(yōu)啟動點的狀態(tài)及軌跡的最大、最小可行啟動高度總結如表2。三個最優(yōu)啟動點的動力著陸軌跡的推力和攻角剖面如圖3所示，圖中時間零點為動力著陸段開始時刻。

圖3 燃料最優(yōu)啟動點的動力著陸軌跡的控制量Fig.3 Control profiles of dynamic landing trajectory of fuel-optimal handover states

表2 燃料最優(yōu)啟動點的狀態(tài)Table 2 The fuel-optimal powered descent initiation states

可見為200 m時，最優(yōu)啟動點的動力著陸軌跡的推力恒為，攻角剖面的形式較為復雜，為100 m和150 m的軌跡中上述性質均成立。為300 m和400 m時，推力剖面為-形式的Bang-Bang控制，即推力由切換至且僅有一次切換，令為推力切換時刻，推力剖面表示為

(12)

為300 m和400 m的軌跡中，攻角剖面近似為Bang-Bang控制，除觸地前較短時間外，攻角取值均為或-，觸地前為滿足終端攻角約束，攻角降至0。其他4條最優(yōu)啟動點不在最小可行高度處的軌跡中該性質均成立，由于此時最優(yōu)啟動點無法通過尋找最小可行啟動高度確定，故下面重點分析為300 m和400 m的軌跡中的性質。

進一步，采用最優(yōu)控制方法研究攻角正負取值與狀態(tài)量的關系。僅考慮縱向平面內的運動，且忽略地球曲率、牽連加速度和Coriolis加速度影響項,和?？紤]攻角較小，引入近似sin≈。將氣動力式(2)代入動力學方程(1)，得動力學方程為

(13)

式(13)對應的哈密爾頓函數為

(14)

推力約束和攻角約束分別如式(5)和式(7)第一式。對上述最優(yōu)控制問題，有定理1如下

設動力著陸問題的最優(yōu)控制為(,)，則必滿足(′-)≤0。

由極小值原理可知(,)使取極小值，將式(14)所示的記為(,)。采用反證法，設(′-)>0，則(,)-(,-)=2(′-)>0，即(,)>(,-)。由于滿足攻角約束式(7)，則-必滿足攻角約束，則(,-)為滿足控制約束的解，則(,)未使取極小值，與“(,)為最優(yōu)控制”矛盾。

數值實驗表明為300 m和400 m的軌跡中最優(yōu)啟動點的動力著陸軌跡中恒有≥0，則由定理1得(′-)≤0。最優(yōu)啟動點不在最小可行高度時，對應的動力著陸軌跡的攻角剖面近似為Bang-Bang控制，僅在觸地前較短時間內攻角降至0。由于攻角由降至0的過程距離觸地的時間較短，此時攻角變化對軌跡影響較小，將動力著陸軌跡中取值均近似為或-，則攻角剖面表示為

(15)

攻角與升力、推力的關系如圖4所示，圖中垂直于速度方向，虛線表示升力，及推力在速度方向的分量，實線表示垂直于方向的力的合力。

圖4 攻角與狀態(tài)量關系Fig.4 Relation between angle of attack and states

本章通過數值優(yōu)化方法分析了返回軌跡中最優(yōu)啟動點的性質。分析表明最優(yōu)啟動點的動力著陸軌跡具有兩種情況，第一種情況為最優(yōu)啟動點的高度為該條軌跡的最小可行高度，此時推力恒為，攻角剖面較為復雜；第二種情況為最優(yōu)啟動高度不在最小可行高度，此時推力剖面如式(12)，攻角剖面近似如式(15)。此時推力剖面中推力切換時刻為待定參數，攻角剖面與狀態(tài)量的關系完全確定。

3 燃料最優(yōu)自適應啟動方法

3.1 基本原理

本文提出的自適應啟動方法通過上述最優(yōu)啟動點的動力著陸軌跡的推力和攻角剖面的性質，在線判斷當前狀態(tài)是否在最優(yōu)啟動點，進而決定是否啟動發(fā)動機。判斷方法為：假設當前狀態(tài)在最優(yōu)啟動點，則存在推力切換時刻，使得以當前狀態(tài)為初始狀態(tài)、推力滿足式(12)，攻角滿足式(15)、側滑角為0的動力著陸軌跡，同時滿足第2章分析中考慮的3個終端約束()=0、()=0、()=。考慮該軌跡的終端時刻自由，即具有兩個未知數和，故松弛1個終端約束，通過2個終端約束確定、，此時動力著陸軌跡唯一確定，若該軌跡滿足松弛的約束，代表當前狀態(tài)在最優(yōu)啟動點。本文松弛約束()=0，通過()=0、()=兩個約束確定、進而確定上述動力著陸軌跡，若該軌跡滿足()近似為0，則啟動發(fā)動機。本文將推力和攻角滿足式(12)和式(15)，滿足約束()=0、()=的軌跡稱為“觸發(fā)軌跡”。

火箭返回中，在氣動減速段的每個制導周期通過觸發(fā)軌跡的()是否接近于0判斷當前狀態(tài)是否在最優(yōu)啟動點。由于氣動減速段火箭自高到低飛行，若最優(yōu)啟動點不在最小可行高度處，且最優(yōu)啟動點的高度小于初始判斷時的高度，則該方法可尋找到最優(yōu)啟動點。本研究發(fā)現，若返回軌跡中水平速度較大，則()隨啟動高度降低而前移。若上一周期()<0且本周期()>0，代表最優(yōu)啟動點在該周期內，故將()>0作為啟動條件，第一次判斷得()>0時啟動發(fā)動機。

本文方法采用的觸發(fā)軌跡的推力和攻角剖面與縱向平面內最優(yōu)啟動點的動力著陸軌跡相同，故在最優(yōu)啟動點處觸發(fā)軌跡的航程與待飛航程相同，避免了文獻[4]和文獻[6]中由于觸發(fā)軌跡與最優(yōu)啟動點的動力著陸軌跡航程存在差異，導致的難以找到最優(yōu)啟動點問題。此外，本方法不依賴于各飛行階段具體的制導方法，僅需要動力著陸段采用燃料最優(yōu)制導。3.2節(jié)將給出上述自適應啟動算法的細節(jié)。

3.2 算法細節(jié)

首先給出觸發(fā)軌跡的計算方法。由于觸發(fā)軌跡的初始狀態(tài)、攻角與狀態(tài)的關系均已給定，觸發(fā)軌跡由推力切換時刻唯一確定，推力切換時刻可通過觸發(fā)軌跡的終端約束()=0、()=求解。令代表滿足上述約束的推力切換時刻，s,代表任意的推力切換時刻。則對任意給定s,，將推力剖面式(12)、攻角剖面式(15)代入動力學方程式(1)進行數值積分，得速度降為時的高度|=。然后，通過Newton-Raphson法迭代求解使得|==0，從而確定觸發(fā)軌跡，并令()代表觸發(fā)軌跡的終端位置。求解s,對應的|=時，由于終端時刻自由，終端條件為=，因此推力切換至后以速度為自變量進行積分。從而將通過()=0、()=求解兩個未知變量,的問題，轉化為通過|==0求解單個未知變量的問題，提高了求解效率。

圖5 整合制導方案及自適應PDI方法流程圖Fig.5 Flowchart of the integrated guidance framework and the proposed adaptive PDI method

4 仿真校驗

4.1 空間質點運動仿真

考慮三組氣動減速段初始狀態(tài)和三個預測制導虛擬目標點，共九條返回軌跡。初始狀態(tài)以表1為標稱值，并考慮初始位置、速度偏差，采用的三組初始狀態(tài)如表3所示。三個虛擬目標點與第2章中分析的三組典型軌跡相同，3自由度模型中三個虛擬目標點的坐標分別為[200, 0, 0]m，[300, 0, 0]m，[400, 0, 0]m，即取上述三個取值，均為0。

表3 氣動減速段初始狀態(tài)Table 3 Initial states of aerodynamic deceleration

9條返回軌跡的燃料消耗與啟動點高度關系如圖6所示，圖中“啟動點”通過提出的自適應啟動方法得到。上述軌跡的啟動點狀態(tài)、最優(yōu)啟動點的燃料消耗和高度，及軌跡的最大、最小可行高度總結如表4，其中啟動時刻以開始氣動減速的時刻為時間零點。可見為200 m時軌跡的最優(yōu)啟動高度為最小可行高度；為300 m和400 m時最優(yōu)啟動高度均不在最小可行高度。從自適應啟動方法確定的啟動點開始動力著陸，上述9條返回軌跡的軌跡剖面如圖7所示。6條含偏差動力著陸軌跡的推力和攻角剖面如圖8所示，3條無偏差軌跡的推力和攻角剖面均與與圖3相近，故圖8中為列出。

圖6 不同初始狀態(tài)下燃料消耗與啟動高度關系圖Fig.6 Profiles of fuel consumption with respect to PDI states under different initial states

圖7 典型返回軌跡的飛行剖面Fig.7 Profiles of the typical return trajectory

圖8 典型返回軌跡啟動點的燃料最優(yōu)動力著陸軌跡的推力和攻角剖面Fig.8 Thethrust and angle of attack profiles of fuel-optimal powered descent trajectories of the triggered states

表4 典型軌跡的啟動點狀態(tài)Table 4 Powered descent initiation states of the typical cases

4.2 計算效率和收斂性

圖9 自適應啟動方法的計算耗時和迭代次數Fig.9 Computational time cost and iterations’ number of theproposed adaptive powered descent initiation method

可見單次求解的最大計算耗時為50.3 ms，最大迭代次數為3次，最大迭代次數出現在首次計算觸發(fā)軌跡時。原因為此時的迭代初值為，隨后計算中由于將上一時刻的作為迭代初值，僅需1至2次迭代即可收斂。

為直觀的展示本文方法的收斂性，做迭代求解時的速度-高度剖面如圖10所示，的迭代初值為-=0。初始狀態(tài)為=300 m時負偏差軌跡的最優(yōu)啟動點的狀態(tài)。圖中表示第步迭代，s,表示該步中切換時刻的取值。圖中第3次迭代即收斂至|==-36×10m，此時滿足迭代收斂條件。

圖10 觸發(fā)軌跡的迭代求解Fig.10 Iterative solutions of trigger trajectories

4.3 蒙特卡洛仿真

為檢驗本文方法關于火箭模型參數的適應性，關于火箭模型參數、初始狀態(tài)進行蒙特卡洛仿真，對比不同參數中自適應啟動方法得到的啟動點與燃料最優(yōu)啟動點間的燃料消耗的差異?？紤]表5所示的9個模型參數和6個初始狀態(tài)參數均在表5所示的取值范圍內服從均勻分布，其中初始縱向位置為相對氣動減速段零攻角、零側滑角航程的偏差，原因為氣動減速段航程受火箭參數和初始狀態(tài)影響較大。火箭的側滑角約束、側滑角變化率約束的取值與攻角相同。氣動減速段制導的虛擬目標點坐標為[300, 0, 0]m，自適應啟動方法的參數與4.1節(jié)相同。

表5 火箭參數和初始狀態(tài)的取值范圍Table 5 Value range of the rocket parameters and initial states

進行1000次蒙特卡洛仿真，自適應啟動方法所得啟動點與GPOPS算法得到的燃料最優(yōu)啟動點間燃料消耗相對差異的統計結果如圖11所示。可見盡管火箭參數和初始狀態(tài)存在較大差異，但本文提出的自適應啟動方法均可得到與燃料最優(yōu)啟動點接近的啟動點。說明盡管本文方法基于特定參數下數值的結論，但在不同火箭參數中具有較好的適應性。

圖11 比最優(yōu)啟動點燃料消耗相對增加量的統計結果Fig.11 Statistical results of the relative fuel usage increase for the triggered PDI states compared with the optimal PDI states

5 結 論

本文針對重復使用火箭垂直返回著陸問題，對動力著陸段的自適應啟動方法展開研究，主要的工作總結如下：

1) 通過數值優(yōu)化方法對返回軌跡系統分析，得到一種解析形式的推力剖面和攻角剖面，描述燃料最優(yōu)啟動點對應的動力著陸軌跡。

2) 設計了一種燃料最優(yōu)的動力著陸段自適應啟動方法。方法基于上述剖面形式判斷當前是否滿足燃料最優(yōu)啟動條件，進而實現自適應啟動。

3) 對提出的自適應啟動方法在多種工況下進行了數值實驗，得到的啟動點與燃料最優(yōu)啟動點十分接近，且迭代求解過程穩(wěn)定，計算效率高。