龍 也，劉一武

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

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

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

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

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

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

0 引 言

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

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

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

1 預(yù)備知識(shí)

1.1 動(dòng)力學(xué)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中，ep0為當(dāng)前能量值，ef為指定的末端能量值.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 仿真分析

3.1 仿真參數(shù)

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

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

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

表3 標(biāo)稱(chēng)進(jìn)入點(diǎn)與理想開(kāi)傘點(diǎn)參數(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è)計(jì)結(jié)果

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

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

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

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

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

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

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

3.3 蒙特卡洛仿真

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

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

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

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

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

4 結(jié) 論

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

［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］龍也，劉一武.火星進(jìn)入段縱向脫敏局限性分析與三自由度脫敏設(shè)計(jì)［J］.宇航學(xué)報(bào)，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—)，男，博士研究生，研究方向?yàn)楹教炱髦茖?dǎo)導(dǎo)航與控制;劉一武(1968—)，男，研究員，研究方向?yàn)楹教炱髦茖?dǎo)導(dǎo)航與控制.

*國(guó)家自然科學(xué)基金資助項(xiàng)目(61403030).

2015-11-29