龍 也，劉一武

(1.北京控制工程研究所，北京100190;2.空間智能控制技術(shù)重點實驗室，北京100190)

一種火星進入段在線脫敏軌跡設(shè)計方法*

龍 也1，2，劉一武1，2

(1.北京控制工程研究所，北京100190;2.空間智能控制技術(shù)重點實驗室，北京100190)

為降低軌跡求解難度，提升脫敏制導(dǎo)的適應(yīng)性，針對火星進入段提出一種在線脫敏軌跡設(shè)計方法.首先，采用預(yù)測的末端航程偏差和狀態(tài)敏感度作為性能指標，利用該指標為傾側(cè)角凸函數(shù)的特性將最優(yōu)求解問題轉(zhuǎn)換為簡單的動態(tài)尋優(yōu)過程;其次，結(jié)合任務(wù)要求和估計的進入點狀態(tài)，通過迭代得到同時滿足航程和橫程要求的三自由度脫敏軌跡.仿真表明該方法可達到與現(xiàn)有脫敏設(shè)計相近的末端狀態(tài)精度.

魯棒制導(dǎo);在線脫敏設(shè)計;火星進入段制導(dǎo);脫敏最優(yōu)控制

0 引 言

目前，火星進入、下降和著陸過程制導(dǎo)、導(dǎo)航與控制技術(shù)的研究聚焦于如何實現(xiàn)0.1 km精確著陸精度［1］.其中，進入段是氣動環(huán)境最惡劣、高度跨度最大，參數(shù)變化及不確定性最多的階段［1］，其末端狀態(tài)精度對實現(xiàn)精確著陸非常關(guān)鍵，研究針對該階段的高精度魯棒制導(dǎo)方法具有重要意義.

進入段制導(dǎo)方法主要分為兩類:標稱軌跡制導(dǎo)和預(yù)測校正制導(dǎo).第一類方法根據(jù)特定性能指標離線優(yōu)化得到標稱軌跡，進入過程中采用反饋線性化+PID［2］、滑模［3］等控制律跟蹤標稱軌跡，對模型不確定性具有一定的魯棒性，但對進入狀態(tài)偏差通常較為敏感.基于脫敏最優(yōu)控制的制導(dǎo)方法［4-6］可較好地解決該問題，它通過在性能指標中添加特定形式的“系統(tǒng)對狀態(tài)擾動敏感度”罰項，獲得具有脫敏特性的標稱軌跡，對進入偏差和過程擾動有較好的魯棒性能，但其敏感度罰項由敏感度傳播方程經(jīng)積分、卷積得到，很難實現(xiàn)軌跡的快速求解.第二類方法采用數(shù)值或解析方法對末端狀態(tài)進行預(yù)測，通過校正過程消除預(yù)測偏差，從而得到可行軌跡和實時控制指令，對進入狀態(tài)偏差魯棒，但迭代求解過程計算量大，需要較高的實時運算能力［7］.

本文在現(xiàn)有脫敏設(shè)計基礎(chǔ)上，綜合標稱制導(dǎo)實時計算量小及預(yù)測校正制導(dǎo)能在線求解可行軌跡的優(yōu)點，針對三自由度脫敏軌跡的在線設(shè)計進行研究.

1 預(yù)備知識

1.1 動力學(xué)模型

忽略自轉(zhuǎn)影響的三自由度火星進入段動力學(xué)方程如下［8］:

式中，r為探測器到火星中心距離，v為速度，γ為航跡傾角，s為待飛航程，θ為經(jīng)度，λ為緯度，ψ為航跡方位角，g=μ/r2為探測器所處位置的重力加速度，μ為火星引力常數(shù)，σ為傾側(cè)角.L和D分別為氣動升力和阻力加速度:

式中，ρ為火星大氣密度，采用標準指數(shù)型大氣模型ρ=ρsexp［－(h/hs)］，h=r－rM為所處位置高度，rM為火星半徑，ρs為火星表面大氣密度，hs為火星大氣尺度高度，CL和CD分別為氣動升力和阻力系數(shù)，Sr為探測器參考面積，m為探測器質(zhì)量.此外，也可采用以下方法計算L、D:

式中，B為彈道系數(shù)，Γ為升阻比系數(shù).

受限于探測器本身物理參數(shù)極限及開傘安全要求，進入段過程中需考慮最大過載、動壓、熱流率等約束［9］:

其中，nmax、qmax、等為前述約束的最大容許值，kq、rn、N、M為熱流率計算相關(guān)系數(shù).

對于標稱制導(dǎo)的軌跡設(shè)計，軌跡除需滿足前述基本過程約束，還需滿足以下基本起終端約束:1)軌跡狀態(tài)初值為標稱進入點狀態(tài);2)軌跡末端狀態(tài)包括高度、速度、航程、經(jīng)緯度(橫縱程)，與預(yù)期值一致.此外，不同制導(dǎo)方法通常采用不同性能指標對滿足基本約束的可行軌跡進行評價，優(yōu)化過程以求解性能最優(yōu)的軌跡為目標.

1.2 現(xiàn)有脫敏設(shè)計

現(xiàn)有脫敏設(shè)計采用縱向或三自由度動力學(xué)方程，但目標函數(shù)中的敏感度罰項均由縱向敏感度傳播方程得出［5-6］:

式中，S(t|t0)=x(t)/x(t0)表示系統(tǒng)狀態(tài)向量x=［rvγs］T的t時刻值對初始值x(t0)的敏感程度，u=cosσ為控制輸入，f(x，u)為式(1)～(4)右側(cè)項式.跟蹤控制律:

式中，u*(t)和x*(t)分別為t時刻的標稱控制量和標稱系統(tǒng)狀態(tài)，K(t)為狀態(tài)反饋增益陣，η(t)為反饋增益乘子.

按式(15)積分得到軌跡各處的敏感度后，根據(jù)不同設(shè)計目標［4-5］計算敏感度罰項，進行離線的軌跡優(yōu)化，即可獲得相應(yīng)脫敏特性的標稱軌跡.

2 在線脫敏軌跡設(shè)計

現(xiàn)有方法的主要問題在于其末端約束的滿足、軌跡敏感度的優(yōu)化均依賴于單次的全局最優(yōu)求解過程.作為對照，文獻［7］采用幅值關(guān)于能量線性的傾側(cè)角剖面對軌跡進行參數(shù)化，求解使末端航程預(yù)測偏差為0的傾側(cè)角，從而得到可行軌跡和實時指令.衍生的阻力加速度制導(dǎo)方法［10］采用三段線性的阻力加速度剖面，通過調(diào)整中間段值滿足末端航程約束，通過迭代航程和反轉(zhuǎn)點得到滿足末端橫程約束的三自由度軌跡.因此，標稱剖面的參數(shù)化以及末端約束的滿足是實現(xiàn)在線設(shè)計的關(guān)鍵.

2.1 預(yù)測方程與軌跡參數(shù)化

預(yù)測方程采用類能量變量e=μ/r－v2/2作為獨立狀態(tài)替代時間t，由縱向動力學(xué)方程式(1)～(4)和敏感度傳播方程式(15)除以=Dv得到:

傾側(cè)角剖面采用文獻［7］方法進行參數(shù)化，由待優(yōu)化的當前傾側(cè)角σp0和固定的末端傾側(cè)角 σf確定:

式中，ep0為當前能量值，ef為指定的末端能量值.

任意軌跡規(guī)劃周期中，系統(tǒng)狀態(tài)已知且σp0確定時，將式(21)代入式(17)～(20)，從ep0積分到ef，即可得到預(yù)測的末端航程spf(σp0)及末端敏感度陣S(tpf|tp0)，tp0為當前時間，tpf為預(yù)測的末端時間.

2.2 性能指標選擇及動態(tài)尋優(yōu)轉(zhuǎn)換

規(guī)劃過程中，算法根據(jù)狀態(tài)約束和性能指標求解出每個規(guī)劃周期的最優(yōu)σp0值.其中，性能指標包含末端敏感度罰項，以優(yōu)化末端對當前狀態(tài)的敏感度，使軌跡具有末端對軌跡各處擾動不敏感特性.狀態(tài)約束僅考慮末端航程，包括過載、熱流率在內(nèi)的其它約束通過更改初始傾側(cè)角和標稱進入狀態(tài)處理［7］.由于傾側(cè)角剖面僅 σp0可調(diào)，航程約束spf(σp0)=0采用罰項形式出現(xiàn)在性能指標中:

式中，λ1≥0為敏感度權(quán)重，J1=|spf|為末端航程偏差罰項，J2為末端敏感度罰項:

式中，Si，j(tpf|tp0)為敏感度矩陣 S(tpf|tp0)的第i行第j列個元素，ci≥0，i=1，2，3，4為權(quán)重因子，表示對各狀態(tài)脫敏特性的注重程度.

圖1 不同能量值下性能指標J與σp0的關(guān)系Fig.1 Relation between performance indexJand σp0corresponding to different energy level

定義歸一化能量E=(e－e0)/(ef－e0)，e0為估計的進入點能量值，由圖1可知，進入過程中，對于E的不同值，J均為σp0的凸函數(shù).因此，規(guī)劃周期fs較小時，結(jié)合傾側(cè)角調(diào)整速率限制，可將逐個周期的最優(yōu)σp0求解問題，轉(zhuǎn)換為如下動態(tài)尋優(yōu)過程:

式中，σp0(N)為上一周期的σp0值，σp0(N+1)為當前σp0尋優(yōu)結(jié)果，為待選性能指標，為最大傾側(cè)角速率.此外，式(24)得到的σp0如為負值則令σp0=0.

由初始傾側(cè)角σ0出發(fā)，按式(24)遞推即可得到縱向軌跡和傾側(cè)角幅值剖面.敏感度罰項的引入會使末端航程出現(xiàn)偏差，即sf≠0，該問題在后續(xù)反轉(zhuǎn)點的迭代求解過程中進行補償.式(24)結(jié)果并非精確的最優(yōu)解，但fs較小時差異可以忽略.

2.3 反轉(zhuǎn)點求解及末端偏差消除

獲得縱向軌跡和傾側(cè)角幅值剖面后，為得到三自由度軌跡，需解決以下問題:1)找到合適的傾側(cè)角反轉(zhuǎn)時間點，消除末端橫向偏差.2)補償敏感度罰項引起的末端航程偏差，并消除縱向偏差.

反轉(zhuǎn)時間點tinv通過牛頓迭代求解，公式為:

式中，tinv(M)為上一次迭代得到的tinv，tinv初始值可設(shè)為tf/2，tf為縱向軌跡的末端時間，為末端橫向偏差，偏導(dǎo)通過有限差分方法近似獲得.進行迭代時，均保留反轉(zhuǎn)點前的傾側(cè)角幅值剖面，按三自由度模型推演至反轉(zhuǎn)點，并按式(24)對后續(xù)軌跡進行三自由度的重規(guī)劃以更新.重規(guī)劃過程中，由式(24)求解得到的σp0均取負，并考慮反轉(zhuǎn)過程中傾側(cè)角的調(diào)整速率限制.

末端航程偏差的補償及縱向偏差的消除，均通過在上述迭代過程中更改初始待飛航程s0實現(xiàn):

3 仿真分析

3.1 仿真參數(shù)

軌跡設(shè)計及蒙特卡洛仿真參數(shù)見表1～3.

表1 相關(guān)物理參數(shù)［11］Tab.1 Related physical parameters

表2 過程約束參數(shù)［1］Tab.2 Process constraint parameters

表3 標稱進入點與理想開傘點參數(shù)［11］Tab.3 Nominal entry point and ideal parachute deploy point parameter

此外K=－［0.01 0.005 50 －0.001］，σ0=0，=0，σf=70°，c1=0.1，c2=0.3，c3=0，c4=1.

表4 蒙特卡洛仿真參數(shù)Tab.4 Monte Carlo simulation parameters

3.2 軌跡設(shè)計結(jié)果

圖2，圖3，表5給出不同λ1值下的縱向軌跡設(shè)計結(jié)果.圖4為λ1=0.8時的三自由度軌跡設(shè)計結(jié)果.表6給出相同計算性能條件下，采用本文與文獻［6］方法分別設(shè)計三自由度軌跡所耗時長的對照結(jié)果.

圖2 不同λ1值時所得縱向軌跡的傾側(cè)角幅值剖面Fig.2 Bank angle magnitude profile of designed longitudinal trajectory corresponding to different λ1values

圖2表明λ1較小時，傾側(cè)角剖面較長時間內(nèi)與僅考慮末端航程時所得剖面基本重合，但一定時間后偏離，λ1增大時偏離時間點提前，λ1=0.95時初期階段即已偏離.圖3表明λ1較小時軌跡末端敏感度改善效果不明顯，λ1過大時末端敏感度反而上升，λ1=0.8時末端敏感度相對較低，說明相對于現(xiàn)有脫敏設(shè)計，本文方法中敏感度權(quán)重并非越大越好.表5表明λ1增大時，軌跡末端高度降低，末端航程偏差增大，因此有必要在設(shè)計過程中進行補償.表6表明本文方法設(shè)計效率相比文獻［6］顯著提升，fs= 0.1 s時僅耗時9分鐘左右，求解更精細軌跡時(fs=0.025 s)耗時半小時左右.

圖3 不同λ1值時所得縱向軌跡的末端敏感度Fig.3 Terminal sensitivity of designed longitudinal trajectory corresponding to different λ1values

表5 不同λ1值時所得縱向軌跡末端狀態(tài)Tab.5 Terminal state of designed longitudinal trajectory corresponding to different λ1values

圖4 λ1=0.8時所得三自由度軌跡的傾側(cè)角幅值剖面Fig.4 Bank angle magnitude profile of designed three degree of freedom trajectory corresponding to λ1=0.8

表6 三自由度軌跡設(shè)計效率對比Tab.6 Efficiency comparison of three degree of freedom trajectory design

3.3 蒙特卡洛仿真

圖5為按文獻［6］方法離線設(shè)計的脫敏軌跡的蒙特卡洛仿真結(jié)果，圖6為λ1=0.8時，本文方法在線設(shè)計的脫敏軌跡的蒙特卡洛仿真結(jié)果.

圖5 離線設(shè)計的三自由度脫敏軌跡蒙特卡洛測試結(jié)果［7］Fig.5 Monte Carlo test results of offline designed 3DOF desensitized trajectory

綜合圖5～6及表7可知，蒙特卡洛測試中，除橫向偏差有一定程度增加外，本文方法與現(xiàn)有脫敏設(shè)計的末端狀態(tài)精度接近.

圖6 在線設(shè)計的三自由度脫敏軌跡蒙特卡洛測試結(jié)果Fig.6 Monte Carlo test results of online designed 3DOF desensitized trajectory

表7 不同方法末端狀態(tài)偏差范圍對照Tab.7 Terminal state deviation scope comparison corresponding to different methods

4 結(jié) 論

現(xiàn)有脫敏設(shè)計需離線求解標稱軌跡，且耗時較長.本文提出一種在線設(shè)計方法，可根據(jù)任務(wù)要求及估計的進入點在線快速設(shè)計三自由度脫敏軌跡.仿真表明該方法具有較高的末端狀態(tài)精度，但所得軌跡末端高度偏低，有待進一步改進.

［1］BRAUN R D，MANNING R M.Mars exploration entry，descent and landing challenges［J］.Journal of Spacecraft and Rockets，2007，44(2):310-323.

［2］TU K Y，MUNIR M S，MEASE K D，Bayard D S.Drag-based predictive tracking guidance for Mars precision landing［J］.Journal of Guidance，Control，and Dynamics，2000，23(4):620-628.

［3］DAI J，XIA Y.Mars atmospheric entry guidance for reference trajectory tracking［J］.Aerospace Science and Technology，2015，45:335-345.

［4］SHEN H J，HANS S，RICHARD W P.Desensitizing the pin-point landing trajectory on Mars［C］//AIAA/ AAS Astrodynamics Specialist Conference and Exhibit.Washington D.C:AIAA，2008.

［5］LI S，PENG Y M.Mars entry trajectory optimization using DOC and DCNLP［J］.Advances in Space Research，2011，47(3):440-452.

［6］龍也，劉一武.火星進入段縱向脫敏局限性分析與三自由度脫敏設(shè)計［J］.宇航學(xué)報，2015，36(8): 861-868.LONG Y，LIU Y W.Limitation analysis of longitudinal desensitize method for Mars entry phase and three degree of freedom desensitize design［J］.Journal of Astronautics，2015，36(8):861-868.

［7］LU P.Predictor-corrector entry guidance for low-lifting vehicles［J］.Journal of Guidance，Control，and Dynamics，2008，31(4):1068-1075.

［8］MANRIQUE J B.Advances in spacecraft atmospheric entry guidance［D］.Irvine:University of California，Irvine，2010.

［9］JOEL B，KENNETH D M.Reachable and controllable sets for planetary entry and landing［J］.Journal of Guidance，Control，and Dynamics，2010，33(3):641-654.

［10］SARAF A，LEAVITT J A，Chen D T，et al.Design and evaluation of an acceleration guidance algorithm for entry［J］.Journal of Spacecraft and Rockets，2004，41(6): 986-996.

［11］LI S，PENG Y M.Command generator tracker based direct model reference adaptive tracking guidance for Mars atmospheric entry［J］.Advances in Space Re-search，2012，49(1):49-63.

An Online Desensitized Trajectory Design Method for Mars Entry Guidance

LONG Ye1，2，LIU Yiwu1，2
(1.Beijing Institute of Control Engineering，Beijing 100190，China;2.Science and Technology on Space Intelligent Control Laboratory，Beijing 100190，China)

To reduce the difficulties in trajectory optimization and improve the adaptability of desensitization guidance，a novel online desensitized trajectory design method is proposed for Mars entry.The weighted value of predicted terminal sail error and state sensitivity are adopted as performance indexes，which are actually convex functions of bank angle.Taking use of this property，the original optimization problem is converted into a simple dynamic optimizing process.Then，according to the mission requirement and estimated entry state，a three-degree-of-freedom desensitized trajectory that simultaneously satisfies the requirement for range and cross range can be found via iteration process.Numerical simulations show that the terminal state accuracy of the proposed method is close to the existing desensitization design.

robust guidance;online desensitization design;mars entry guidance;desensitized optimal control

V448.2

:A

:1674-1579(2016)02-0020-06

10.3969/j.issn.1674-1579.2016.02.004

龍 也(1990—)，男，博士研究生，研究方向為航天器制導(dǎo)導(dǎo)航與控制;劉一武(1968—)，男，研究員，研究方向為航天器制導(dǎo)導(dǎo)航與控制.

*國家自然科學(xué)基金資助項目(61403030).

2015-11-29