權(quán)申明 , 王竹 晁濤 楊明

(1. 哈爾濱工業(yè)大學(xué) 控制與仿真中心，黑龍江 哈爾濱150080; 2. 上海機(jī)電工程研究所，上海 201109)

0 引言

飛行器在再入點(diǎn)的狀態(tài)如再入速度傾角、再入位置、再入速度等，對(duì)飛行器在大氣層內(nèi)的飛行狀態(tài)和飛行能力有至關(guān)重要的影響[1-2]。良好的再入狀態(tài)，對(duì)飛行器在大氣層內(nèi)的駐點(diǎn)熱流密度、總加熱量等物理約束、滑翔飛行的調(diào)節(jié)能力產(chǎn)生積極的作用，同時(shí)能夠增加飛行器的安全性和穩(wěn)定性，減輕制導(dǎo)控制系統(tǒng)的壓力。再入點(diǎn)各狀態(tài)約束值由飛行器設(shè)計(jì)的總體單位根據(jù)任務(wù)需求及飛行器性能計(jì)算得到，設(shè)計(jì)離軌段制導(dǎo)控制算法時(shí)可認(rèn)為是已知的。

飛行器離軌制動(dòng)段制導(dǎo)算法可以規(guī)劃出一條滿足再入點(diǎn)狀態(tài)的轉(zhuǎn)移軌道。為方便研究，在初步設(shè)計(jì)軌道時(shí)，通常假設(shè)飛行器在軌道轉(zhuǎn)移過(guò)程中可瞬間獲得所需的速度沖量且僅受地球中心引力場(chǎng)作用。

文獻(xiàn)[3]指出，良好的離軌制動(dòng)飛行狀態(tài)是確保天地往返飛行器安全、準(zhǔn)確返回既定著陸場(chǎng)區(qū)的前提。傳統(tǒng)鈍頭體再入返回艙，受到橫向機(jī)動(dòng)能力和搜救能力限制，必須保證再入點(diǎn)有足夠的控制精度。離軌制動(dòng)可以視為軌道轉(zhuǎn)移的一種形式，現(xiàn)有的軌道轉(zhuǎn)移軌跡規(guī)劃算法設(shè)計(jì)時(shí)，發(fā)動(dòng)機(jī)主要有脈沖推力(視為瞬時(shí)脈沖)、小推力和有限推力三種形式。其中，基于脈沖推力的計(jì)算較為簡(jiǎn)單，相對(duì)成熟的Lambert 制導(dǎo)算法可在較短時(shí)間內(nèi)計(jì)算出變軌點(diǎn)的速度增量，然而實(shí)際飛行過(guò)程中，發(fā)動(dòng)機(jī)難以提供任意大小的發(fā)動(dòng)機(jī)推力；小推力軌道轉(zhuǎn)移通常用在深空探測(cè)等場(chǎng)景，推力較小，因此需要數(shù)日才能完成某任務(wù)；有限推力離軌制動(dòng)形式在航天器返回、載荷釋放等場(chǎng)景下應(yīng)用廣泛，對(duì)該問(wèn)題的研究具有工程意義。具體地，有限推力離軌規(guī)劃又可分為逼近式策略以及優(yōu)化式策略。

在逼近式策略方面，文獻(xiàn)[4]提出兩次“推-滑”的離軌制導(dǎo)策略，分析了推力沿不同方向施加時(shí)能量和動(dòng)量矩的變化規(guī)律，給出臨界高度的計(jì)算公式及兩次“推-滑”的判斷條件。文獻(xiàn)[5]將離軌制動(dòng)過(guò)程分成能量耗散段和閉路導(dǎo)引控制段，探討了在燃料隨機(jī)耗盡情況下推力方向?qū)υ偃朦c(diǎn)參數(shù)的影響，分析了推力方向切換與能量窗口的關(guān)系。文獻(xiàn)[6]利用軌道瞬時(shí)根數(shù)進(jìn)行離軌制動(dòng)算法設(shè)計(jì)，將離軌制動(dòng)問(wèn)題轉(zhuǎn)化為制動(dòng)起點(diǎn)的緯度幅角的單參數(shù)搜索問(wèn)題，并考慮了J2 項(xiàng)攝動(dòng)、初始狀態(tài)偏差、質(zhì)量與推力誤差的影響。文獻(xiàn)[7]借鑒運(yùn)載火箭能量管理的思路，為實(shí)現(xiàn)耗盡關(guān)機(jī)，將離軌制動(dòng)過(guò)程分為能量耗散與閉環(huán)制導(dǎo)控制兩個(gè)階段，在運(yùn)算效率與制導(dǎo)精度方面均有一定的提升。逼近式策略方面的研究基于機(jī)理分析，從能量角度設(shè)計(jì)制動(dòng)方案，通常僅考慮再入點(diǎn)速度大小和方向約束，由于沒(méi)有考慮時(shí)間代價(jià)，傳統(tǒng)的逼近式策略難以滿足快速離軌需求。

在優(yōu)化式策略方面，文獻(xiàn)[8]提出的方法能夠?qū)崿F(xiàn)指定再入點(diǎn)位置的飛行器離軌任務(wù)，但在內(nèi)層計(jì)算過(guò)程中需要反復(fù)調(diào)用求解非線性規(guī)劃求解器，因此難以實(shí)現(xiàn)制動(dòng)速度的在線計(jì)算。文獻(xiàn)[9]考慮再入返回飛行器對(duì)再入點(diǎn)高度和再入角的要求，推算了離軌制動(dòng)的軌道參數(shù)及所需速度增量，同時(shí)指出了慣性系速度傾角和地固系下的再入角差異較小。 文獻(xiàn)[10]針對(duì)建立了有限推力轉(zhuǎn)移軌跡優(yōu)化的凸優(yōu)化模型，并提出一種迭代逼近算法，但每步迭代計(jì)算中的凸優(yōu)化求解耗時(shí)在1.24~3.23 s 之間，算法難以實(shí)現(xiàn)在線優(yōu)化。文獻(xiàn)[11-12]針對(duì)燃料最優(yōu)問(wèn)題進(jìn)行研究，在飛行器總體優(yōu)化設(shè)計(jì)階段具有重要意義，能夠節(jié)約大量成本。文獻(xiàn)[13]基于最優(yōu)控制理論推導(dǎo)了離軌機(jī)動(dòng)的最優(yōu)性條件，給出了協(xié)態(tài)變量和狀態(tài)變量的解析解，進(jìn)而轉(zhuǎn)化為兩點(diǎn)邊值問(wèn)題，但初值的選取對(duì)求解結(jié)果有較大影響。文獻(xiàn)[14]針對(duì)轉(zhuǎn)移時(shí)間最優(yōu)和燃料最優(yōu)兩種性能指標(biāo)，將原問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題，將Bezier 曲線法得到的結(jié)果作為初值，采用內(nèi)點(diǎn)法與序列二次規(guī)劃算法進(jìn)行求解，雖然對(duì)初值的計(jì)算進(jìn)行了改進(jìn)與優(yōu)化，但是仍難以做到在線計(jì)算。文獻(xiàn)[15]分別針對(duì)共面、異面軌道轉(zhuǎn)移問(wèn)題，研究了基于凸優(yōu)化方法的軌跡規(guī)劃，并選用不同指標(biāo)進(jìn)行了軌道機(jī)動(dòng)能力分析，但是算法實(shí)時(shí)性較差，難以工程實(shí)現(xiàn)。文獻(xiàn)[16]提出一種組合制導(dǎo)方法，前期采用燃料最優(yōu)的開環(huán)策略，采用基于標(biāo)稱軌跡的閉環(huán)制導(dǎo)，受開環(huán)控制的影響，該方法的調(diào)節(jié)能力有待進(jìn)一步研究。文獻(xiàn)[17]針對(duì)動(dòng)能彈垂直再入問(wèn)題利用自適應(yīng)偽譜法對(duì)有限推力制動(dòng)與連續(xù)小推力制動(dòng)兩種方式進(jìn)行了研究。上述研究成果往往需要借助工具箱，如SNOPT、SQP、CVX 等[18]。

由以上分析可以看出，基于逼近式和優(yōu)化式策略的離軌制動(dòng)算法均存在一定的不足：逼近式策略算法簡(jiǎn)單，計(jì)算效率高，易于工程實(shí)現(xiàn)，然而少有考慮時(shí)間最優(yōu)的離軌需求；優(yōu)化式策略可以求解同時(shí)考慮多個(gè)再入狀態(tài)約束下的快速離軌制動(dòng)問(wèn)題，然而優(yōu)化計(jì)算效率較低，難以實(shí)時(shí)在線計(jì)算。隨著空間飛行器智能化需求的不斷提高，再入飛行任務(wù)已經(jīng)從到達(dá)地面固定目標(biāo)點(diǎn)，逐步向著覆蓋全球任意位置發(fā)展，因此傳統(tǒng)的固定點(diǎn)離軌制動(dòng)、固定的再入難以適應(yīng)這種需求。對(duì)于已經(jīng)部署就位的航天器，其上所攜帶燃料固定，為了應(yīng)對(duì)一定范圍內(nèi)可變的任務(wù)要求，飛行器需要具備一定的再入點(diǎn)狀態(tài)調(diào)節(jié)能力。

本文研究在離軌制動(dòng)任務(wù)確定之后，燃料有限情況下的多終端約束的時(shí)間最優(yōu)離軌制動(dòng)算法。本文的主要?jiǎng)?chuàng)新點(diǎn)如下：

1)結(jié)合優(yōu)化算法的結(jié)果，通過(guò)其最優(yōu)軌道形式，提出燃料約束情況下的有限推力快速離軌制動(dòng)制導(dǎo)算法；

2)提出虛擬再入角的概念，將時(shí)間最優(yōu)快速離軌制動(dòng)問(wèn)題參數(shù)化為虛擬再入角的求解，簡(jiǎn)化優(yōu)化問(wèn)題；

3)為滿足離軌制導(dǎo)算法的實(shí)時(shí)性要求，本文所提出的離軌制動(dòng)算法計(jì)算簡(jiǎn)單、快速且具有較強(qiáng)的收斂性。

1 問(wèn)題描述

飛行器軌道轉(zhuǎn)移示意圖如圖1所示。圖1中，M為飛行器接到任務(wù)時(shí)的位置，K為離軌制動(dòng)初始位置，E為再人點(diǎn)位置，C為地球表面目標(biāo)點(diǎn)，h為軌道高度，hf為離軌段終端高度。選取軌道轉(zhuǎn)移時(shí)間固定或能量最小等為指標(biāo)后，可唯一確定出轉(zhuǎn)移軌道。

圖1 軌道轉(zhuǎn)移示意圖Fig. 1 Schematic diagram of orbital transfer

由于推力有限且連續(xù)作用，飛行狀態(tài)不能根據(jù)自由段解析彈道的特征解析計(jì)算，而是需要通過(guò)積分彈道方程來(lái)獲得。飛行器在大氣層外的運(yùn)動(dòng)模型可以描述為

式中：r為飛行器位置矢量，r為飛行器質(zhì)心距地心的徑向距離，r=‖r‖；v為飛行器速度矢量；μ為地球引力常數(shù)；pe為發(fā)動(dòng)機(jī)推力，pe表示制動(dòng)推力大小；m為飛行器質(zhì)量；Isp為發(fā)動(dòng)機(jī)比沖；g0為地球表面重力加速度大小。

為了更加清晰地描述飛行器相對(duì)于地球的運(yùn)動(dòng)情況，給出慣性系下飛行器運(yùn)動(dòng)模型的另一種表示：

式中：v為飛行器相對(duì)地速；θ為飛行器速度傾角，在上為正；λ和φ分別為飛行器星下點(diǎn)的經(jīng)度和緯度；χ為彈道偏角，順時(shí)針為正；α表示制動(dòng)角，即推力方向與速度方向的夾角(類似于攻角的定義)；g為重力加速度，g=g0(R0/r)2，R0為地球平均半徑。

2 離軌制動(dòng)算法設(shè)計(jì)

優(yōu)化式策略通?；谧顑?yōu)控制相關(guān)理論，求解復(fù)雜、效率低，往往需要調(diào)用優(yōu)化工具箱得到最優(yōu)結(jié)果；逼近式策略基于能量守恒等原理，結(jié)合衛(wèi)星軌道理論，實(shí)時(shí)計(jì)算可行的控制量，逐漸調(diào)整軌跡，可逐步逼近期望終端狀態(tài)。

2.1 沖量假設(shè)的離軌制動(dòng)算法

沖量假設(shè)簡(jiǎn)化了飛行器離軌制動(dòng)的運(yùn)動(dòng)過(guò)程，認(rèn)為飛行器可瞬間獲得需要速度，給離軌制動(dòng)段的初步設(shè)計(jì)帶來(lái)了方便。

針對(duì)更加一般性離軌制動(dòng)問(wèn)題，飛行器一旦進(jìn)入離軌制動(dòng)階段，為實(shí)現(xiàn)高度的降低，飛行軌跡為橢圓軌道。

設(shè)v1為當(dāng)前位置需用速度，r1表示由地心指向飛行器實(shí)時(shí)位置的向量，r2表示由地心指向期望再入點(diǎn)位置的向量，h0為飛行器離軌點(diǎn)高度，θ1為根據(jù)能量守恒和動(dòng)量矩不變?cè)矸辞蟪鰎1處的期望速度傾角。沖量假設(shè)下的離軌制動(dòng)問(wèn)題描述為：已知r1處的初始速度v0和初始速度傾角θ0，r2處的終端速度v2和終端速度傾角θf(wàn)，求解沖量假設(shè)下的離軌段制動(dòng)參數(shù)Δθ和Δv，v1和制動(dòng)速度的夾角為φ。矢量r1和r2的大小分別為r1和r2，v1和v2的大小分別為v1和v2。圖2 為橢圓軌道離軌段航跡示意圖。

圖2 橢圓軌道離軌段航跡示意圖Fig. 2 Deorbit trajectory of an elliptical orbit

下面針對(duì)圖2 所示的任意橢圓軌道進(jìn)行算法設(shè)計(jì)。根據(jù)能量守恒和動(dòng)量矩不變?cè)?，?/p>

通過(guò)式(3)可以求解出θ1和v1，代入下式中：

由此可以得到離軌段制動(dòng)參數(shù)Δθ和Δv，若能夠瞬時(shí)沿著φ的方向提供一個(gè)Δv，則飛行器將沿著新軌道運(yùn)動(dòng)，最終到達(dá)再入點(diǎn)。

2.2 有限推力離軌制動(dòng)算法

通常情況下，天基再入滑翔飛行器的離軌制動(dòng)發(fā)動(dòng)機(jī)推力有限，因此2.1 節(jié)的沖量假設(shè)離軌制動(dòng)制導(dǎo)方法難以適用。此外，即使發(fā)動(dòng)機(jī)可瞬時(shí)提供任意方向及大小的推力，由于空間中攝動(dòng)因素的存在，其離軌段終端將產(chǎn)生較大偏差，影響再入飛行。因此，采用有限推力方式的離軌制動(dòng)制導(dǎo)方法研究更有實(shí)際意義，脈沖假設(shè)下的計(jì)算結(jié)果可作為參考。

由于燃料的消耗，飛行器質(zhì)量不斷減少，根據(jù)給定的發(fā)動(dòng)機(jī)參數(shù)，可以利用推力和比沖進(jìn)行計(jì)算。

在有限推力模型中，采用沖量假設(shè)下的離軌制動(dòng)制導(dǎo)方法，可以求解得到?jīng)_量假設(shè)下的離軌段離軌參數(shù)Δθ和Δv。在求解有限推力下的離軌段軌跡時(shí)，給定天基再入滑翔飛行器的制動(dòng)點(diǎn)狀態(tài)及再入點(diǎn)狀態(tài)約束，計(jì)算推力大小與方向變化曲線，控制飛行器按照期望再入角、再入速度到達(dá)再入點(diǎn)高度。

結(jié)合2.1 節(jié)瞬時(shí)脈沖仿真結(jié)果，設(shè)計(jì)一種采用逼近式策略的考慮有限推力的制導(dǎo)步驟如下：

步驟1根據(jù)飛行器實(shí)時(shí)狀態(tài)、再入角度及速度，由式(3)按照瞬時(shí)脈沖方式確定出當(dāng)前時(shí)刻所需的速度和角度信息。

步驟2結(jié)合飛行器當(dāng)前角度和速度信息，計(jì)算制動(dòng)速度Δv及制動(dòng)角度。

步驟3沿著制動(dòng)角度方向，改變飛行器速度增量為Δv，發(fā)動(dòng)機(jī)采用最大推力。

步驟4若所需速度小于其設(shè)定閾值，即，則關(guān)閉發(fā)動(dòng)機(jī)，飛行器沿橢圓軌道自由運(yùn)動(dòng)，否則循環(huán)步驟1~步驟3 的計(jì)算。

2.3 快速離軌制動(dòng)算法設(shè)計(jì)

通過(guò)分析不難發(fā)現(xiàn)，以上兩種方法沒(méi)有考慮時(shí)間最優(yōu)性指標(biāo)，僅考慮了再入點(diǎn)高度、速度大小與速度傾角，因此求得的離軌制動(dòng)軌跡難以滿足快速離軌的需求。

通過(guò)最優(yōu)控制理論中的時(shí)間最優(yōu)結(jié)論可知，時(shí)間最優(yōu)情況下控制量通常是Bang-Bang 控制，類似于開關(guān)控制量，離軌制動(dòng)階段飛行器所攜帶的燃料有限，因此假設(shè)時(shí)間最優(yōu)的離軌制動(dòng)問(wèn)題中，制動(dòng)發(fā)動(dòng)機(jī)推力大小與制動(dòng)角度形式如圖3 所示。圖3中，O點(diǎn)為原點(diǎn)，α軸表示制動(dòng)角，t軸為時(shí)間，α1(t)為0~t1時(shí)間內(nèi)的制動(dòng)角變化曲線，α2(t)為t2~t3時(shí)間內(nèi)的制動(dòng)角變化曲線，pemax為最大制動(dòng)推力，t1為制動(dòng)發(fā)動(dòng)機(jī)首次關(guān)機(jī)時(shí)刻，t2為第2 次開機(jī)時(shí)刻，t3為第2 次關(guān)機(jī)時(shí)刻，t4為到達(dá)再入點(diǎn)的時(shí)刻。

圖3 時(shí)間最優(yōu)制動(dòng)發(fā)動(dòng)機(jī)工作狀態(tài)示意圖Fig. 3 Schematic diagram of time-optimal deorbit

求解α(t)，使得終端時(shí)間tf最小，同時(shí)滿足初始各狀態(tài)與終端各狀態(tài)，可描述為

式中：u為制動(dòng)角變化曲線α(t)；Δm為飛行器總質(zhì)量變化量；mfuel為燃料質(zhì)量。

可將式(5)描述的問(wèn)題轉(zhuǎn)化為非線性規(guī)劃問(wèn)題，進(jìn)而基于偽譜法進(jìn)行求解，但是該類方法計(jì)算效率較差，工程實(shí)現(xiàn)性不強(qiáng)。還可將該非線性規(guī)劃問(wèn)題轉(zhuǎn)化為凸優(yōu)化問(wèn)題，凸優(yōu)化在求解最優(yōu)控制問(wèn)題的快速性與最優(yōu)性上均有顯著提升[19]，然而現(xiàn)階段凸優(yōu)化算法仍難以滿足在線實(shí)時(shí)計(jì)算的速度需求。

由于傳統(tǒng)的有限推力離軌制動(dòng)算法并未考慮時(shí)間最優(yōu)性，存在一定不足。設(shè)計(jì)改進(jìn)的有限推力離軌制動(dòng)算法如圖 4 虛線曲線所示，引入虛擬再入角的概念，在前期采用虛擬再入角，后期采用實(shí)際期望再入角作為約束。

圖4 改進(jìn)方法與傳統(tǒng)方法彈道傾角變化曲線示意圖Fig. 4 Curve of the flight path angle with the improved method and the traditional method

定義1虛擬再入角：

在再入返回的初期，設(shè)計(jì)不同于實(shí)際期望終端速度傾角θf(wàn)的某固定或者時(shí)變?cè)偃虢亲鳛橹茖?dǎo)輸入指令，稱角度為虛擬再入角(在后期仍使用θf(wàn)作為期望再入角)。

在離軌制動(dòng)前期主要控制再入速度傾角與再入速度大小，后期根據(jù)文獻(xiàn)[4]中得到的垂直速度方向施加力不改變能量大小的結(jié)論，通過(guò)預(yù)測(cè)剩余燃料消耗時(shí)間，進(jìn)行正推與反推的能量耗散策略。后期使用剩余燃料朝著目標(biāo)點(diǎn)飛行，由于后期時(shí)間相對(duì)短，終端速度不會(huì)偏差太大。本文設(shè)計(jì)的快速離軌制動(dòng)制導(dǎo)仿真流程如圖 5 所示。

圖5 快速離軌制動(dòng)制導(dǎo)仿真流程圖Fig. 5 Flow chart of simulation of fast deorbit guidance

因此，可將問(wèn)題式(5)轉(zhuǎn)化為如下問(wèn)題：

并使用凸優(yōu)化方法對(duì)關(guān)鍵參數(shù) Δt1、Δt2，以及第一階段的虛擬期望終端傾角求解。

基于沖量假設(shè)的有限推力下離軌制動(dòng)算法存在的問(wèn)題在于沒(méi)有考慮快速離軌的需求，但在再入速度傾角和速度大小的控制上，計(jì)算量小且精度較高。基于凸優(yōu)化的離軌制動(dòng)可以考慮時(shí)間最優(yōu)約束，仍難以實(shí)時(shí)計(jì)算，但是通過(guò)與有限推力離軌制動(dòng)算法相結(jié)合，將最優(yōu)結(jié)果作為有限推力制動(dòng)方法虛擬終端約束，從而給出快速離軌制動(dòng)制導(dǎo)算法。

2.4 基于凸優(yōu)化方法的子問(wèn)題求解

文獻(xiàn)[19]給出了一種可以處理非線性等式的凸優(yōu)化問(wèn)題的方法。在時(shí)間最優(yōu)的離軌制動(dòng)問(wèn)題中，給定起始時(shí)刻的高度、速度，要求飛行器以最短時(shí)間到達(dá)某一高度，并滿足終端的速度大小及傾角 要求。

使用凸優(yōu)化對(duì)問(wèn)題進(jìn)行求解需要分為三個(gè)步驟。

2.4.1 狀態(tài)方程歸一化

式中：σ為松弛變量，代表推力-加速度。設(shè)狀態(tài)量為x=(r,v,z)T，控制量為u=(τ,σ)T，將狀態(tài)方程寫為如下矩陣形式：

當(dāng)使用迭代方法求解最優(yōu)解時(shí)，式(10)可看作線性方程。

2.4.2 求解初始軌跡

使用迭代方法求解時(shí)間最優(yōu)離軌制動(dòng)問(wèn)題時(shí)，首先需要一條初始軌跡，在此基礎(chǔ)進(jìn)行優(yōu)化。

求解時(shí)間最優(yōu)離軌問(wèn)題，需要將離軌過(guò)程分為N段，每一段時(shí)間間隔為t0/N，求解N+ 1個(gè)狀態(tài)變量及控制變量。飛行器初始位置、速度與終端時(shí)刻位置、速度約束可表示為

離軌過(guò)程約束為式(9)~式(13)，由于地心距r為歸一化后的變量，在離軌過(guò)程中近似為1，在對(duì)結(jié)果精確度要求不高的情況下，式(9)中的r可近似為1，變?yōu)榫€性方程。

由于所有約束及目標(biāo)函數(shù)均為凸的，此問(wèn)題為一標(biāo)準(zhǔn)的凸優(yōu)化問(wèn)題。

2.4.3 迭代求解時(shí)間最優(yōu)離軌制動(dòng)問(wèn)題

求解時(shí)間最優(yōu)離軌制動(dòng)問(wèn)題不指定固定的再入點(diǎn)以及離軌時(shí)間，只要求到達(dá)固定的再入高度時(shí)滿足指定的速度大小及傾角，因此需要引入時(shí)間變量。此時(shí)式(9)與式(13)均為非線性等式約束。同時(shí)，為保證結(jié)果的精確，需要在初始解的基礎(chǔ)上進(jìn)行多次迭代求解。

當(dāng)凸優(yōu)化問(wèn)題中含有非線性等式約束hi(y) = 0,i= 1，… ,q時(shí)，假設(shè)當(dāng)前迭代結(jié)果為y[k]，y= （f x,u,tf）。求解y[k+1]的過(guò)程如下：

在y[k]處對(duì)函數(shù)進(jìn)行泰勒展開，hi(y) =0可以寫為

同時(shí)，又有展開式

凸優(yōu)化問(wèn)題的約束為

優(yōu)化結(jié)果為yp，將在y[k]處1階展開：

式(18)兩側(cè)同乘(yp-y[k])T/2，得到

因此便將非凸約束hi(y) =0轉(zhuǎn)換為凸約束。

時(shí)間最優(yōu)的離軌制動(dòng)問(wèn)題中，除了考慮式(9)~式(11)的約束之外，其余約束為

以離軌段終端時(shí)間tf最短為指標(biāo)。式(9)的離散形式以及式(22)為非線性約束，將其按前述可得到凸函數(shù)的約束條件，并附加約束

即可求得y[k+1]以及u[k+1]，多次迭代直至收斂，得到時(shí)間最優(yōu)離軌制動(dòng)結(jié)果。

至此，將離軌制動(dòng)問(wèn)題轉(zhuǎn)化為凸優(yōu)化問(wèn)題，使用2階錐規(guī)劃方法進(jìn)行求解，求解工具詳見文獻(xiàn)[20]。

2.5 基于最優(yōu)軌跡的虛擬指令計(jì)算

根據(jù)凸優(yōu)化求解得到的優(yōu)軌跡中發(fā)動(dòng)機(jī)推力大小曲線，可以看出時(shí)間最優(yōu)的離軌制動(dòng)問(wèn)題中，控制量呈現(xiàn)“推-滑-推”的特點(diǎn)，從而根據(jù)其時(shí)間變化規(guī)律可得 Δt1、 Δt2，根據(jù)速度、傾角變化曲線，根據(jù)式(4)可以反求出虛擬再入角θ?f大小。

3 數(shù)值仿真及分析

考慮有兩臺(tái)推力發(fā)動(dòng)機(jī)，每臺(tái)推力為500 N，比沖為300 s，離軌到再入階段單獨(dú)計(jì)算燃料質(zhì)量，燃料為250 kg。

根據(jù)2.1~2.3 節(jié)給出的不同方法，考慮從某圓軌道進(jìn)行離軌制動(dòng)，初始高度h0=500 km，初始速度為v0=7 616.7 m/s，初始傾角θ0= 0°，期望終端高度hf=120 km、終端速度vf=7 950 m/s、終端傾角θf(wàn)=-2 .5°為例，分別給出不同方法的仿真結(jié)果。需要指出的是，本文所提出的快速離軌制動(dòng)制導(dǎo)方法采用類似的策略在橫向進(jìn)行虛擬指令設(shè)計(jì)后，可適用于三維空間中的異面離軌制動(dòng)，仿真中僅以縱向運(yùn)動(dòng)為例。

3.1 基于凸優(yōu)化仿真結(jié)果

將快速離軌制動(dòng)問(wèn)題轉(zhuǎn)化為凸優(yōu)化問(wèn)題，進(jìn)而使用2階錐規(guī)劃進(jìn)行求解。仿真結(jié)果如圖6所示。通過(guò)多約束下時(shí)間最優(yōu)離軌制動(dòng)的仿真計(jì)算，可以看出飛行軌跡能夠?qū)崿F(xiàn)時(shí)間最優(yōu)的指標(biāo)；從推力大小曲線可以看出，發(fā)動(dòng)機(jī)工作狀態(tài)類似Bang-Bang控制，與最優(yōu)控制的時(shí)間最優(yōu)控制問(wèn)題結(jié)果一致；由圖6(f)可以看出，飛行器終端時(shí)刻燃料用盡，終端各狀態(tài)滿足要求。同時(shí)，為了驗(yàn)證方法有效性，放寬燃料限制，當(dāng)燃料足夠時(shí)，推力可以長(zhǎng)時(shí)間最大推力工作，通過(guò)改變制動(dòng)角來(lái)實(shí)現(xiàn)時(shí)間最優(yōu)指標(biāo)。

圖6 多約束下時(shí)間最優(yōu)離軌制動(dòng)狀態(tài)量變化曲線Fig. 6 Curves of time-optimal deorbit with multiple constraints

綜合以上仿真曲線可以看出，時(shí)間最優(yōu)下的仿真結(jié)果雖比有限推力離軌制動(dòng)的制動(dòng)時(shí)間指標(biāo)更優(yōu)，但是從編程實(shí)現(xiàn)及運(yùn)算效率上來(lái)看，代價(jià)較大。即使在初值合理的情況下，單次優(yōu)化耗時(shí)在 1 s 以上，因此難以實(shí)現(xiàn)在線制導(dǎo)的目標(biāo)。

3.2 快速離軌制動(dòng)算法仿真

為驗(yàn)證本文算法的有效性，選取不同終端狀態(tài)θf(wàn)=-1 .5° 、θf(wàn)=-2 .5° 、θf(wàn)=-3 .5° 、θf(wàn)=- 4.5° 進(jìn)行仿真，基于文獻(xiàn)[6]中在線制導(dǎo)算法的仿真結(jié)果如圖 7 所示，基于本文快速離軌制動(dòng)制導(dǎo)算法仿真結(jié)果如圖 8 所示。

對(duì)比圖7、圖8仿真結(jié)果曲線可以看出，兩種方法均能滿足終端高度、速度、傾角的終端約束。在傳統(tǒng)算法中，未考慮快速離軌制動(dòng)需求時(shí)，燃料有剩余，且隨著終端傾角(幅值)的增大，離軌制動(dòng)段飛行時(shí)間減小，燃料消耗增加，這是由于空間軌道中，速度方向的改變需要較大的速度增量。對(duì)比以上圖7(f)、圖8(f)子圖燃料質(zhì)量變化曲線可知，為滿足離軌制動(dòng)快速性的需求，本文算法能夠?qū)⑷剂虾谋M，離軌段節(jié)約時(shí)間最大值約為500 s。在計(jì)算效率方面，單步計(jì)算耗時(shí)均在5 ms以內(nèi)，能夠滿足在線制導(dǎo)的需求。

圖7 傳統(tǒng)離軌制動(dòng)制導(dǎo)算法仿真結(jié)果Fig. 7 Simulation results of the traditional deorbit guidance algorithm

圖8 快速離軌制動(dòng)制導(dǎo)算法仿真結(jié)果Fig. 8 Simulation results of the fast deorbit guidance algorithm

3.3 蒙特卡洛仿真

在實(shí)際離軌制動(dòng)過(guò)程中，飛行器發(fā)動(dòng)機(jī)受到各種偏差因素的影響，本文主要考慮推力大小攝動(dòng)Δpe、推力方向攝動(dòng)Δφ、燃料質(zhì)量攝動(dòng)Δm，各偏差均為均勻分布如表1 所示。

表1 參數(shù)攝動(dòng)表Table 1 Parameter perturbations

進(jìn)行1 000 組蒙特卡洛仿真實(shí)驗(yàn)，仿真結(jié)果如圖 9 所示(Rf表示離軌制動(dòng)段射程)。

圖9 蒙特卡洛仿真參數(shù)散布圖Fig. 9 Distribution of Monte Carlo simulation parameters

根據(jù)終端狀態(tài)散布結(jié)果可以看出，該算法對(duì)終端速度與終端傾角控制精度較高，其中終端速度偏差為±0.2 m/s，終端傾角偏差為±0.002°。由于該算法未對(duì)射程進(jìn)行精確控制，該狀態(tài)的絕對(duì)偏差較大為±50 km，然而該值相對(duì)于整個(gè)再入階段相對(duì)誤差較小，約為0.5%，同時(shí)考慮到后續(xù)的再入滑翔階段，對(duì)射程調(diào)節(jié)能力較強(qiáng)，因此離軌制動(dòng)階段射程的影響可忽略。

4 結(jié)論

本文針對(duì)快速離軌制動(dòng)的需求，提出基于虛擬再入角的快速離軌制動(dòng)制導(dǎo)方法。借鑒固體運(yùn)載火箭助推段末期的能量管理策略，通過(guò)設(shè)置合理的虛擬再入角與虛擬指令控制邏輯，設(shè)計(jì)“做正功”與“做負(fù)功”的切換規(guī)律，進(jìn)行多余能量的耗散，可以實(shí)現(xiàn)快速離軌制動(dòng)算法的在線設(shè)計(jì)。得出主要結(jié)論如下：

1)引入虛擬再入角的概念將離軌制動(dòng)最優(yōu)控制問(wèn)題參數(shù)化，采用少次的凸優(yōu)化方法求解最優(yōu)軌跡進(jìn)行虛擬再入角的更新，該算法相比于傳統(tǒng)凸優(yōu)化方法具有計(jì)算簡(jiǎn)單、運(yùn)算效率高的特點(diǎn)。

2)相比于現(xiàn)有制導(dǎo)算法，基于虛擬再入角的快速離軌制動(dòng)算法考慮了離軌耗時(shí)的因素，在燃料耗盡的情況下，能夠節(jié)約最多500 s 的離軌時(shí)長(zhǎng)，滿足了離軌制動(dòng)快速性的需求。

3)該算法對(duì)離軌階段發(fā)動(dòng)機(jī)推力大小與方向、燃料質(zhì)量攝動(dòng)具有較強(qiáng)的魯棒性，再入點(diǎn)速度控制精度為±0.2 m/s，再入角控制精度為±0.002°。

在本文工作的基礎(chǔ)上，后續(xù)將考慮J2 攝動(dòng)的快速離軌制動(dòng)在線制導(dǎo)算法，考慮地球橢球體引起的非球形引力的等因素的影響，進(jìn)一步增強(qiáng)本文算法的工程意義。