周晚萌，李海陽，王 華，黃海兵，李云飛

(1.國防科技大學空天科學學院應用力學系，長沙410073；2.長沙翔宇信息科技有限公司，長沙410073；3.上海宇航系統(tǒng)工程研究所，上海201100)

1 引言

近年來，隨著國際空間站服役期限的臨近，原空間站合作方開始共同提出了一項關(guān)于國際環(huán)月空間站的研究計劃，截止于2016年10月，研究方案已基本明確。同時，我國也在開展圍繞環(huán)月空間站的載人探月系列任務，以充分繼承我國在近地空間站建設的經(jīng)驗優(yōu)勢，實現(xiàn)月球的長期探測。載人登月系統(tǒng)復雜，很難由某家單位獨立實現(xiàn)全任務仿真，必須利用仿真平臺進行多家任務單位的模型集成。載人登月軟著陸任務仿真系統(tǒng)主要包括飛行任務仿真單元、飛行任務規(guī)劃單元、飛行任務評估單元和可視化單元。

飛行任務仿真單元對月面著陸器和上升飛行器組合體在下降和返回過程進行動力學仿真。再入動力學模型一般采用速度、航跡角以及偏航角描述速度狀態(tài)，但月面軟著陸的最終速度為0，求解過程中會出現(xiàn)奇異[1]，為此飛行任務仿真單元直接利用速度分量來描述動力學方程中的速度狀態(tài)。飛行任務規(guī)劃單元主要針對定點著陸和月面上升問題展開研究。這類問題均可以歸結(jié)為終端時間可變的燃料最優(yōu)控制問題，目前主要分為直接法和間接法兩類[2]。間接法將通過極大值原理建立協(xié)調(diào)方程與狀態(tài)方程聯(lián)立求解[3-4]，但該方法對協(xié)調(diào)變量初值要求較高[5]，同時軟著陸過程中的約束也加大了方程的推導難度[6]，因此國內(nèi)外主要采用直接法對此類問題展開研究。直接法具有更強的魯棒性和適應性[7]，直接法中的直接打靶法將控制量離散，并將問題轉(zhuǎn)化為多設計參數(shù)的非線性規(guī)劃問題，利用SQP算法或智能優(yōu)化算法來求解[8-10]，主要問題是直接法為了降低求解難度，忽略了月球自轉(zhuǎn)，因此并未引入著陸點的位置約束。另外，也有學者利用直接法中的正交配點法求解月球軟著陸問題[11-13]。

美國、俄羅斯/蘇聯(lián)、歐洲和日本多年來致力于仿真平臺的構(gòu)建，并在地面試驗驗證中起到了非常重要的作用[14-15]。我國很多高校和研究機構(gòu)在相關(guān)領域也搭建了各具特色的仿真平臺，如吳魁等[16]采用分布式仿真架構(gòu)提出了運載火箭總體設計仿真評估平臺；孫福煜等[17]設計了基于高層體系架構(gòu)(High Level Architecture，HLA)的分布仿真環(huán)境KD-HLA；卿杜政等[18]提出了基于HLA的SSS-RTI系統(tǒng)。

本文采用共享內(nèi)存與HLA相結(jié)合的機制與可視化組裝技術(shù)，搭建了從建模、調(diào)試、組裝到仿真的全生命周期的數(shù)字化載人登月任務仿真平臺(Simulation Management Platform，SIM)。 在 SIM平臺的基礎上建立由飛行任務規(guī)劃單元、飛行任務仿真單元、飛行任務評估單元和可視化單元組成的載人登月軟著陸任務仿真系統(tǒng)。利用該系統(tǒng)進行軟著陸飛行任務規(guī)劃算法的驗證和仿真評估，其中飛行任務規(guī)劃單元采用變網(wǎng)格Radau偽譜法實現(xiàn)，以期在應急條件約束[19]下，可以同時保證計算精度和收斂效率。

2 飛行任務仿真單元

2.1 動力學模型

如圖1所示，飛行仿真單元采用月面著陸坐標系o-xyz建立動力學模型，著陸坐標系屬于動坐標系，具體定義為z軸沿飛行器矢徑方向，xoy平面為當?shù)厮矫妫瑇軸與當?shù)卣狈较虼嬖谝粋€偏角γ，且保持不變，y軸由右手定則確定。推力器的推力方向可以由方位角β和高低角α來確定。高低角定義為推力矢量與當?shù)厮矫娴膴A角，范圍為 [-90°，90°]， 正值表示在水平面上方，負值表示在水平面下方。方位角定義為推力水平面投影矢量與著陸坐標系的夾角，范圍是 [0 °，360°]。 為獲得完整的動力學模型，首先需要確定月心J2000坐標系O-XYZ、月心赤道固聯(lián)坐標系O-X1Y1Z1與月面著陸坐標系三者之間的轉(zhuǎn)換關(guān)系。僅考慮月球自轉(zhuǎn)影響，O-XYZ與O-X1Y1Z1之間相對角速度為月球自轉(zhuǎn)角速度ωL，O-X1Y1Z1與o-xyz的角速度為Ω，需要先繞Z1軸轉(zhuǎn)動一個經(jīng)度角λ，再繞新的Y1軸轉(zhuǎn)動90°-φ，φ為緯度角，最后繞新的Z1軸轉(zhuǎn)動180°-γ。

圖1 月面著陸坐標系示意圖Fig.1 Coordinate systems in lunar soft-landing

月球自轉(zhuǎn)角速度ωL、o-xyz相對于O-X1Y1Z1的角速度為Ω，動力學模型的o-xyz分量可以表示為式(1)[20]。

其中，質(zhì)心相對于O-X1Y1Z1的速度在o-xyz坐標系的分量為(u，v，w)，r為著陸器矢徑，T為推力幅值，m為質(zhì)量，μL為月球引力常數(shù)。另外，全過程燃料消耗外加質(zhì)量方程為式(2)。

其中，Isp表示發(fā)動機比沖，g0表示地球重力加速度。為提高軌跡優(yōu)化精度，將距離單位DU，時間單位TU，質(zhì)量單位MU等單位量級歸一化，具體定義 1 DU＝1RL，1 TU ＝，1 MU＝1m0。 其中，RL為月球半徑，m0為組合體初始質(zhì)量。

2.2 邊界條件及過程約束

軟著陸的邊界條件包括位置和速度約束，一種位置約束為式(3)。

其中，h為距月面高度，i＝0、f表示初始和終端狀態(tài)。速度的初始約束表示為式(4)。

其中，V0表示初始軌道的速度，由于組合體動力下降之前處于遠月點高度為環(huán)月軌道高度hLLO、近月點高度為h0的橢圓軌道，并從近月點開始降軌，終端約束為u(tf)＝0，v(tf)＝0，w(tf)＝0。

在動力下降過程中還要考慮到其他的過程約束，發(fā)動機推力以及轉(zhuǎn)角的限制為式(5)。

航天員承受過載限制為式(6)。

其中，過載n可以用最大值的情況來直接表示為n＝T/m＋μL/r2。

要滿足應急返回的條件，必須使軟著陸過程中的速度分量以及徑向與切向速度比值不能超過限制，如式(7)所示。

3 飛行任務規(guī)劃單元

按照配點的選取方式不同，可以將偽譜法分為Legendre、Gauss以及Radau偽譜法，這里采用Radau偽譜法作為基礎，并進一步動態(tài)劃分網(wǎng)格。

Radau配點為-1，1( ]上多項式gK(τ)的根，稱為Legendre-Gauss-Radau(LGR)點，多項式表示為式(8)。

其中，PK(τ)為K階勒讓德多項式，可表示為式(9)。

將時間變量歸一化為τ-∈ [-1，1]，以τ-＝-1和K個LGR點作為插值節(jié)點，對狀態(tài)變量進行K階拉格朗日插值，如式(10)所示。

其中，Li(τ-)為拉格朗日插值基函數(shù)，表示為式(11)。

由式(10)可知，插值函數(shù)在節(jié)點處的真實值等于近似值，即(i＝0，…K)。同理，可得控制變量的插值結(jié)果如式(12)所示。

其中，為拉格朗日插值基函數(shù)，表示為式(13)。

離散化后在LGR點上的狀態(tài)微分，可以表示為式(14)。

其中，k＝1，…，K，D∈R RK×(K＋1)為微分矩陣，與優(yōu)化問題本身無關(guān)，只與LGR點的選取有關(guān)，為拉格朗日插值基函數(shù)的微分，表達式為式(15)。

其中，Gk(τ)＝(1＋τ)gk(τ)。 由上述的狀態(tài)微分近似，在LGR點(配點)上對動力學微分方程進行配置，如式(16)所示。

系統(tǒng)的終端狀態(tài)可表示為式(17)。

最后，對性能指標J＝Φ(X0，t0，Xf，tf)＋進行處理，對其中的拉格朗日型指標采用高斯求積進行離散，可得式(19)。

偽譜法計算過程中，配點并非全局插值，而是將求解區(qū)間劃分為多個小區(qū)間，每個區(qū)間內(nèi)具有多個插值，一般的偽譜法在每個區(qū)間內(nèi)以相同的階數(shù)對狀態(tài)和控制變量進行插值，而變網(wǎng)格自適應則要求通過式(16)計算子區(qū)間殘差， 判斷是否需要進一步加密網(wǎng)格或調(diào)整插值階數(shù)。若在容許誤差范圍內(nèi)，則停止迭代計算，否則需要自適應調(diào)整。網(wǎng)格細化的判據(jù)為式(20)。

其中，為第k個區(qū)間內(nèi)的最大曲率和平均曲率，區(qū)間內(nèi)第s個點的曲率計算公式為式(21)。

若r(k)超過給定的最大限度，則需要按式(22)加密網(wǎng)格；否則，需按式(23)提高插件階數(shù)。

式中，A、B為調(diào)整參數(shù)。

在每次迭代中，網(wǎng)格的劃分和插值階的選取都不斷趨于合理，在保證計算精度和速度的條件下，增強了對不連續(xù)和狀態(tài)量快變問題的優(yōu)化計算能力。

4 仿真系統(tǒng)設計

在硬件平臺的基礎上，構(gòu)建應用層、工具層、中間件層、基礎資源庫層4層結(jié)構(gòu)的軟件環(huán)境。應用層針對特定任務(如載人登月軟著陸任務)開發(fā)仿真應用，工具層、中間件層、基礎資源庫層和硬件層構(gòu)成實驗室通用的仿真運行環(huán)境，如圖2所示。工具層中的仿真運行支撐平臺SIM作為載人航天任務預研仿真論證的主要平臺。

4.1 SIM平臺架構(gòu)設計

圖2 SIM平臺仿真運行環(huán)境架構(gòu)Fig.2 Structure of simulation environment in SIM platform

SIM平臺包括仿真管理軟件SIM、仿真客戶端軟件SIM.Client、仿真數(shù)據(jù)管理軟件SIM.Data、仿真模型生成軟件SIM.Creator、仿真模型調(diào)試軟件SIM.Debug等，如圖3所示。為了解決傳統(tǒng)HLA仿真完全依賴網(wǎng)絡交互數(shù)據(jù)，難以適應高頻度數(shù)據(jù)交互要求的難題，SIM仿真平臺采用高性能HLA/SharedMemory兩層交互機制，在同一計算機上的模型之間采用SharedMemory交互機制，在不同計算機上的模型之間采用HLA交互機制，極大提高了仿真速度。

在SIM平臺中，支持仿真模型框架生成(SIM.Creator)、仿真模型調(diào)試(SIM.Debug)、仿真模型集中管理(SIM/SIM.Client)以及數(shù)據(jù)管理(SIM.Data)。

平臺采用“分布計算、集中管理”思想，在管理節(jié)點集中進行仿真模型管理和仿真試驗管理。仿真試驗管理形成仿真想定后，由仿真分布管理工具根據(jù)仿真想定將各仿真模型和模型參數(shù)分布至相應計算節(jié)點，并將模型啟動、關(guān)閉等操作集中在管理節(jié)點進行，如圖4所示。

圖3 SIM平臺組成Fig.3 Composition of SIM platform

4.2 任務仿真系統(tǒng)設計

載人登月軟著陸任務仿真系統(tǒng)采用SIM作為底層平臺，整個系統(tǒng)分為4個單元：飛行任務規(guī)劃單元、飛行任務仿真單元、任務評估單元和可視化單元。

圖4 仿真分布管理模型分發(fā)運行原理圖Fig.4 Schematic diagram of model operation and dispersal in simulation distribution management

在對載人登月軟著陸任務仿真之前，需通過飛行任務規(guī)劃單元按照飛行任務要求，對載人登月軟著陸涉及的軟著陸軌跡進行規(guī)劃，并將初始軌道、變軌控制參數(shù)等規(guī)劃結(jié)果作為初始條件的一部分發(fā)送至飛行任務仿真單元。通過SIM.Creator自動生成月面著陸器飛行仿真單元、任務評估單元以及可視化單元的模型框架，建立仿真模型。利用SIM.Debug獨立于仿真平臺進行單節(jié)點調(diào)試并錄入到SIM。SIM采用可視化的樹狀結(jié)構(gòu)管理各節(jié)點的錄入模型，通過拖拽、連續(xù)等方式實現(xiàn)試驗系統(tǒng)的圖形化組裝。集中配置試驗參數(shù)和仿真模型，完成仿真系統(tǒng)的設計。SIM平臺真正實現(xiàn)了貫穿仿真系統(tǒng)的全生命周期過程，在建模、調(diào)試運行、仿真評估及態(tài)勢展現(xiàn)的各個階段提供全方位支持。根據(jù)飛行任務規(guī)劃結(jié)果，完成飛行任務仿真推演，實現(xiàn)軟著陸全過程的仿真模擬，實時將飛行器與相關(guān)設備的狀態(tài)與參數(shù)輸出至飛行任務評估單元和可視化單元，完成測控通信、光照能源等專項評估以及可視化展示，如圖5所示。

5 仿真分析

將搭建好的軟著陸飛行任務規(guī)劃單元和飛行任務仿真單元分布部署在不同的仿真節(jié)點上。節(jié)點采用的處理器為3.4 GHz Intel Core i7-6700，8GBRAM，操作系統(tǒng)為Windows7。

采用算例的具體方案為：組合體飛行器從環(huán)月極軌道出發(fā)，進入高度為15 km×110 km的橢圓軌道，無動力滑行至近月點后，開始進行動力下降。

圖5 載人登月軟著陸任務仿真系統(tǒng)架構(gòu)圖Fig.5 Architecture of manned lunar soft-landing simulation system

軟著陸過程中的推力最小值Tmin為4.5 kN，最大值Tmax為45 kN。 高低角范圍-30°～30°，方位角范圍160°～200°。 最大過載nmax為1，vmax為50 m/s，wmax為45 m/s，Kmax為0.3。

變網(wǎng)格Radau偽譜法的各子區(qū)間的配點個數(shù)取值范圍是4～10，迭代的最大代數(shù)為20，容許誤差為10-5。按照表1的數(shù)據(jù)配置各單元模型，開始運行仿真，同時記錄Gauss偽譜法和變網(wǎng)格Radau偽譜法最終的燃料消耗與計算耗時，結(jié)果如表2所示。

表1 仿真算例參數(shù)配置Table 1 Parameters of simulation example

表2 仿真結(jié)果對比Table 2 Comparison of simulation results

由此，可知本文算法規(guī)劃結(jié)果燃料為5953 kg，與Gauss偽譜法相比，燃料消耗更少。對比計算耗時可以發(fā)現(xiàn)變網(wǎng)格Radau偽譜法僅為98.32 s，低于Gauss偽譜法的平均計算耗時。由于變網(wǎng)格的劃分以及插值階數(shù)的改變僅發(fā)生在局部網(wǎng)格而非全局過程，因此所提方法在保證計算結(jié)果的同時還可以獲得高于傳統(tǒng)Gauss偽譜法的計算效率。飛行過程中的推力方向如圖6所示。

圖6 推力方向曲線Fig.6 Directions of the thrust

由圖中的仿真結(jié)果可知，推力的方位角為鈍角，并且位于水平面上方，始終保持縱向和法向的減速。在中間一段平穩(wěn)的飛行段，推力方向約為180°并持續(xù)了200 s左右的時間。軟著陸過程中的推力幅值如圖7所示。

圖7 推力幅值與過載大小Fig.7 Thrust amplitude vs.overloads

圖7 中左邊y軸展示了推力幅值的大小，右邊y軸展示了整個過程的過載大小。從圖中可知，全過程推力器始終處于開機狀態(tài)，推力恒定為最大值，過載在約束范圍內(nèi)不斷升高。軟著陸過程中的剩余質(zhì)量如圖8所示，速度分量如圖9所示。

圖8 燃料消耗曲線Fig.8 Fuel consumption

圖9 組合體速度分量Fig.9 Velocity components of the spacecraft complex

從圖中可以看出，全程處于開機狀態(tài)，剩余質(zhì)量在不斷減少，同時3個速度分量最終達到終端速度要求。切向速度不斷減小，而橫向和徑向由于存在應急條件約束，在中間時段速度大小保持在邊界值附近，持續(xù)時間約為200 s。全過程的哈密爾頓函數(shù)結(jié)果如圖10所示。整個軟著陸過程的軌跡如圖11所示。

圖10 哈密爾頓函數(shù)結(jié)果Fig.10 Results of Hamiltonian function

圖11 組合體軟著陸飛行軌跡Fig.11 Lunar soft-landing trajectory

由圖10可知，全過程的哈密爾頓函數(shù)值約為0，滿足一階最優(yōu)的必要條件，偽譜法能保證計算結(jié)果的最優(yōu)性。從軌跡結(jié)果可知，整個軟著陸過程軌跡平滑，無過大的機動，具備實際工程價值。

6 結(jié)論

1)基于SIM平臺搭建的載人登月軟著陸任務仿真系統(tǒng)采用月面著陸坐標系下的速度分量建立任務仿真單元，同時采用變網(wǎng)格Radau偽譜法作為任務規(guī)劃單元，能夠很好地實現(xiàn)軟著陸任務過程的仿真。

2)變網(wǎng)格Radau偽譜法與其他方法的計算結(jié)果對比可知，算法可以在滿足計算精度的同時，提高計算效率并獲得了更優(yōu)的計算結(jié)果。

3)由仿真結(jié)果可知，軟著陸過程中總速度不斷減少，且推力矢量始終位于水平面上方。在軟著陸過程中橫向和法向的速度幅值會先達到最大限制并以該極限值飛行一段時間后再逐漸降為0。