解永鋒，周文勇，杜大程，陳 益，王傳魁

（北京宇航系統(tǒng)工程研究所，北京，100076）

有限推力機(jī)動(dòng)變軌能力包絡(luò)快速計(jì)算分析

解永鋒，周文勇，杜大程，陳 益，王傳魁

（北京宇航系統(tǒng)工程研究所，北京，100076）

通過設(shè)計(jì)加權(quán)的組合性能指標(biāo)函數(shù)，將機(jī)動(dòng)變軌能力包絡(luò)計(jì)算評(píng)估問題轉(zhuǎn)化為一系列使組合性能指標(biāo)最優(yōu)的飛行任務(wù)規(guī)劃問題求解，并采用求解精度高、收斂速度快的偽譜結(jié)點(diǎn)-非線性規(guī)劃方法將存在多個(gè)推進(jìn)-滑行段的最優(yōu)任務(wù)規(guī)劃問題離散轉(zhuǎn)換并快速求解。算例仿真表明，采用該方法可以綜合評(píng)估限定推進(jìn)劑耗量、程序角速率限幅等約束下的機(jī)動(dòng)變軌能力包絡(luò)，給出滿足可控性要求的最優(yōu)飛行軌道，進(jìn)而可以增強(qiáng)小推力長(zhǎng)期在軌機(jī)動(dòng)飛行器的多樣化任務(wù)適應(yīng)能力。

機(jī)動(dòng)變軌；能力包絡(luò)；任務(wù)規(guī)劃；偽譜結(jié)點(diǎn)

0 引 言

為滿足空間多樣化任務(wù)需求，軌道機(jī)動(dòng)飛行器在軌機(jī)動(dòng)過程通常分為變軌任務(wù)規(guī)劃和主動(dòng)變軌實(shí)施兩個(gè)階段，在變軌任務(wù)規(guī)劃階段，需要根據(jù)飛行器當(dāng)前飛行軌道以及目標(biāo)軌道進(jìn)行變軌策略規(guī)劃，包括對(duì)發(fā)動(dòng)機(jī)工作次數(shù)、點(diǎn)火時(shí)刻、主動(dòng)段飛行程序角規(guī)劃；在主動(dòng)變軌實(shí)施階段，需要根據(jù)規(guī)劃方案，進(jìn)行主動(dòng)變軌前調(diào)姿，然后軌控發(fā)動(dòng)機(jī)按要求工作，實(shí)現(xiàn)變軌機(jī)動(dòng)。對(duì)于長(zhǎng)期在軌采用有限推力的軌道機(jī)動(dòng)飛行器而言，由于攜帶的推進(jìn)劑有限，必須在綜合評(píng)估其剩余機(jī)動(dòng)變軌能力的基礎(chǔ)上，根據(jù)任務(wù)目標(biāo)，制定兼顧機(jī)動(dòng)變軌時(shí)間和燃料耗量的變軌策略，以保障軌道機(jī)動(dòng)飛行器的后續(xù)任務(wù)執(zhí)行能力。

目前有關(guān)空間軌道轉(zhuǎn)移優(yōu)化的文獻(xiàn)較多[1～4]，但大多是針對(duì)固定目標(biāo)軌道的優(yōu)化算法研究，通常采用簡(jiǎn)化的軌道動(dòng)力學(xué)模型，不考慮程序角速率限幅和程序角限制，以燃料消耗最優(yōu)作為性能指標(biāo)函數(shù)，優(yōu)化得到的飛行軌跡并未考慮工程可實(shí)現(xiàn)性，而優(yōu)化算法亦不能用于真實(shí)飛行任務(wù)規(guī)劃和變軌能力評(píng)估層面。本文以最優(yōu)控制理論為基礎(chǔ)，采用偽譜結(jié)點(diǎn)-非線性規(guī)劃方法，通過引入偽控制量，將機(jī)動(dòng)變軌能力包絡(luò)計(jì)算評(píng)估問題轉(zhuǎn)化為一系列使組合性能指標(biāo)最優(yōu)的機(jī)動(dòng)變軌任務(wù)規(guī)劃問題，規(guī)劃得到的飛行軌道程序角平滑且可自動(dòng)滿足角速率限幅約束，具備工程可實(shí)現(xiàn)性。

1 空間軌道動(dòng)力學(xué)描述

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

地心發(fā)射慣性坐標(biāo)系下建立的空間飛行器的三自由度質(zhì)點(diǎn)運(yùn)動(dòng)方程為

式中r,V為慣性位置和速度向量；g(r)為考慮 J2引力項(xiàng)攝動(dòng)的重力加速度；T為發(fā)動(dòng)機(jī)推力，與飛行器縱軸方向相同；vacT為最大真空推力幅值；η為發(fā)動(dòng)機(jī)節(jié)流閥值，η＞0；spI為發(fā)動(dòng)機(jī)比沖；m為飛行器質(zhì)量；g0為海平面重力加速度； 1b為單位向量，表示縱軸方向，與俯仰及偏航程序角(,)φ ψ的對(duì)應(yīng)關(guān)系為

1.2 飛行約束

機(jī)動(dòng)變軌需滿足的狀態(tài)變量和控制量約束條件為

式中minmax,ηη為節(jié)流閥約束的上、下邊界，若發(fā)動(dòng)機(jī)為固定推力，ηmin=ηmax=1；φmin,φmax,ψmin,ψmax為俯仰、偏航程序角約束的上、下邊界；,為程序角速率

限幅值。

2 偽譜最優(yōu)控制

2.1 最優(yōu)控制問題描述

針對(duì)空間軌道規(guī)劃問題的特點(diǎn)，定義最優(yōu)控制問題(B)為：a）尋求狀態(tài)-控制函數(shù)對(duì)(x, u)∈RNx×RNu和可能未知的端點(diǎn)時(shí)刻τ0和τf；b）以最小化設(shè)定的性能指標(biāo)函數(shù)；c）滿足狀態(tài)-控制變量約束、動(dòng)態(tài)約束（狀態(tài)方程約束）、軌道約束和端點(diǎn)約束等，表示為

并且滿足：

式中x為狀態(tài)變量，x∈RNx；xL,xU分別為狀態(tài)變量約束的上、下界；u為控制變量，u∈RNu；uL,uU分別為控制變量約束的上、下界；E為終端性能指標(biāo)，；F，f為性能指標(biāo)函數(shù)，函數(shù)中的拉格朗日項(xiàng)、狀態(tài)方程右端項(xiàng)為，；h為狀態(tài)-控制混合軌道約束，；hL,hU分別為軌道約束的上、下邊界；e為端點(diǎn)約束，；uL,uU分別為端點(diǎn)約束的上、下邊界。

2.2 偽譜最優(yōu)控制問題轉(zhuǎn)化

為達(dá)到節(jié)省燃料的目的，機(jī)動(dòng)變軌可能以推進(jìn)-滑行-推進(jìn)的模式進(jìn)行，對(duì)于這類多階段并存在狀態(tài)量、控制量不連續(xù)的最優(yōu)控制問題，可采用偽譜結(jié)點(diǎn)法[5]進(jìn)行處理，將最優(yōu)控制問題轉(zhuǎn)換為非線性規(guī)劃問題，并利用成熟的非線性規(guī)劃方法求解。

不失一般性，以兩階段最優(yōu)控制問題為例，每一階段的問題形式與2.1節(jié)問題()B相同。定義τe是中間時(shí)間點(diǎn)（偽譜結(jié)點(diǎn)），兩段時(shí)間子區(qū)間分別為，則最優(yōu)控制問題可以離散轉(zhuǎn)化為

以及事件時(shí)間點(diǎn) τ0,τe,τf，以滿足如下約束：

式中N1, N2為計(jì)算節(jié)點(diǎn)數(shù)；為(Ni+1)×(Ni+1)維微分矩陣，為N階拉格朗日插值多項(xiàng)式，并最小化性能指標(biāo)為

關(guān)于偽譜結(jié)點(diǎn)法的詳細(xì)介紹可參考文獻(xiàn)[5]，在此不再詳述。

2.3 光滑控制生成

為保證控制量-俯仰及偏航程序角變化平滑，引入程序角的變化率作為偽控制量，而原控制量作為新的狀態(tài)變量，擴(kuò)展后的狀態(tài)方程為

2.4 性能指標(biāo)

若要衡量變軌機(jī)動(dòng)對(duì)軌道面的改變能力，則可定義軌道傾角及升交點(diǎn)赤經(jīng)的加權(quán)組合作為性能指標(biāo)函數(shù)即：

式中w為權(quán)重因子，w∈[0,1]；if,?f分別為主動(dòng)變軌段終端軌道傾角及升交點(diǎn)赤經(jīng)。并限定發(fā)動(dòng)機(jī)最大工作時(shí)間以表征所允許的最大推進(jìn)劑耗量。軌道傾角及升交點(diǎn)赤經(jīng)同時(shí)調(diào)整的能力包絡(luò)計(jì)算可以轉(zhuǎn)化為求解給定終端飛行約束和路徑約束，并采用考慮J2引力項(xiàng)攝動(dòng)的軌道動(dòng)力學(xué)方程，以最優(yōu)化上述組合性能指標(biāo)的最優(yōu)控制問題。

對(duì)于終端軌道固定的機(jī)動(dòng)變軌任務(wù)，為使燃料消耗最優(yōu)，可采用如下性能指標(biāo)函數(shù)：

3 算例分析

3.1 異面機(jī)動(dòng)軌道規(guī)劃

假設(shè)在軌機(jī)動(dòng)飛行器初始質(zhì)量1 860 kg、主發(fā)動(dòng)機(jī)推力2 000 N、比沖為290.5 s，推進(jìn)劑加注量700 kg，飛行器初始及終端軌道約束見表1。俯仰及偏航程序角限幅值取1.5（°）/s。以燃料消耗最優(yōu)作為性能指標(biāo)，分別基于偽譜最優(yōu)控制方法規(guī)劃一次變軌和兩次變軌模式的飛行軌道。

表1 機(jī)動(dòng)變軌段端點(diǎn)約束

3.1.1 一次變軌模式

圖1 俯仰程序角及其變化率

圖2 偏航程序角及其變化率

飛行軌道參數(shù)變化曲線如圖3所示。

圖3 軌道傾角及升交點(diǎn)赤經(jīng)

由圖3可以看出，一次變軌同時(shí)調(diào)整俯仰及偏航程序角，可以達(dá)到所要求的多終端約束，進(jìn)入目標(biāo)橢圓軌道。通過對(duì)圖 3的分析可知，由于所需修正的升交點(diǎn)赤經(jīng)較大且軌道傾角修正量很小，而主發(fā)動(dòng)機(jī)點(diǎn)火時(shí)刻軌道緯度幅角約為40.82°，由軌道攝動(dòng)理論可知，對(duì)升交點(diǎn)赤經(jīng)修正的同時(shí)必然會(huì)不可避免地造成軌道傾角的改變，為達(dá)到所要求的目標(biāo)修正量，必須在提升近地點(diǎn)高度的同時(shí)，拉大當(dāng)前軌道傾角與目標(biāo)值的偏差，然后再實(shí)現(xiàn)軌道傾角與赤經(jīng)的同時(shí)修正，這也是造成推進(jìn)劑耗量較大的主要原因。

3.1.2 兩次變軌模式

采用前面所述偽譜結(jié)點(diǎn)法，將兩次主動(dòng)段工作時(shí)間及滑行時(shí)間作為優(yōu)化變量，在兩次主動(dòng)段之間加入滑行段，其中第1主動(dòng)段主要用于改變半長(zhǎng)軸和偏心率，將軌道拉起，然后經(jīng)過一段時(shí)間的滑行，待緯度幅角滿足要求后，發(fā)動(dòng)機(jī)再次工作，對(duì)軌道面進(jìn)行修正。采用兩次變軌模式，可顯著減小推進(jìn)劑耗量，主發(fā)動(dòng)機(jī)兩次工作時(shí)間之和為340.5 s，推進(jìn)劑消耗239.1 kg。

兩次變軌段的飛行程序角見圖4、圖5，與一次變軌模式相比，程序角變化范圍較小，更易于控制實(shí)現(xiàn)。

圖4 第1主動(dòng)段飛行程序角

圖5 第2主動(dòng)段飛行程序角

與之對(duì)應(yīng)的軌道傾角及升交點(diǎn)赤經(jīng)變化見圖6、圖7，與一次變軌模式相比，軌道傾角變化幅角較小，且升交點(diǎn)赤經(jīng)逐步向目標(biāo)值靠近，不會(huì)出現(xiàn)先增大再減小的情況。

圖6 軌道傾角

圖7 升交點(diǎn)赤經(jīng)

3.2 軌道修正能力包絡(luò)

采用2.4節(jié)定義的性能指標(biāo)函數(shù)，采用的初始及終端約束見表 2，計(jì)算經(jīng)過一次機(jī)動(dòng)變軌進(jìn)入155 km× 200 km 的橢圓軌道，并同時(shí)調(diào)整軌道傾角及升交點(diǎn)赤經(jīng)能力包絡(luò)。

表2 機(jī)動(dòng)變軌段端點(diǎn)約束

因推進(jìn)劑加注量為700 kg，考慮主發(fā)動(dòng)機(jī)最大允許工作時(shí)間tbmax在200～800 s之間變化，計(jì)算得到的異面變軌機(jī)動(dòng)能力包絡(luò)曲線如圖8所示。不同線形表示不同的允許主發(fā)動(dòng)機(jī)工作時(shí)間，所組成的閉合區(qū)域表示可以達(dá)到的軌道機(jī)動(dòng)量，且所需主發(fā)動(dòng)機(jī)工作時(shí)間不超過tbmax，由圖8可見所能實(shí)現(xiàn)的軌道可機(jī)動(dòng)區(qū)域隨 tbmax的減小而縮?。辉谲壍纼A角改變量|iΔ|較小的情況下，主發(fā)動(dòng)機(jī)工作時(shí)間（推進(jìn)劑耗量）隨升交點(diǎn)赤經(jīng)改變量|?Δ|的增大而快速增大。圖中A點(diǎn)表示表 1所列的主發(fā)動(dòng)機(jī)點(diǎn)火時(shí)刻目標(biāo)軌道與飛行器的軌道傾角及升交點(diǎn)赤經(jīng)之差，靠近tbmax=600 s閉合區(qū)域邊界線，在tbmax=700 s閉合區(qū)域內(nèi)，與優(yōu)化計(jì)算得到的變軌時(shí)間相吻合。

圖8 異面變軌能力包絡(luò)

4 結(jié) 論

采用基于偽譜最優(yōu)控制的軌道規(guī)劃策略，可以快速評(píng)估限定推進(jìn)劑耗量、程序角速率限幅等多重約束下的機(jī)動(dòng)變軌能力包絡(luò)，并給出滿足工程可控性要求的一次或多次變軌最優(yōu)飛行軌道，進(jìn)而增強(qiáng)有限推力長(zhǎng)期在軌軌道機(jī)動(dòng)飛行器的多樣化任務(wù)適應(yīng)能力。

[1]Yue X, Yang Y, Geng Z. Indirect optimization for finite-thrust time-optimal orbital maneuver[J]. Journal of Guidance Control and Dynamics, 2010, 33(2): 628-634.

[2]王常虹, 曲耀斌, 陸智俊, 等. 航天器有限推力軌道轉(zhuǎn)移的軌跡優(yōu)化方法[J]. 西南交通大學(xué)學(xué)報(bào), 2013, 48(2): 390-394.

[3]鄭昌, 文趙偉. 基于進(jìn)化計(jì)算的軌道轉(zhuǎn)移時(shí)間-能量?jī)?yōu)化方法[J]. 計(jì)算機(jī)工程與設(shè)計(jì), 2008, 29(3): 735-737.

[4]劉峰. 有限推力軌道轉(zhuǎn)移優(yōu)化方法與應(yīng)用研究[D]. 長(zhǎng)沙: 國(guó)防科技大學(xué), 2009.

[5]Ross I M, Fahroo F. Discrete verification of necessary conditions for switched nonlinear optimal control systems[C]. Boston: Massachusetts, 2004.

Rapid Computation and Analysis of Finite Thrust Maneuver Capability

Xie Yong-feng, Zhou Wen-yong, Du Da-cheng, Chen Yi, Wang Chuan-kui
(Beijing Institute of Astronautical Systems Engineering, Beijing, 100076)

Evaluating of maneuver capability envelope is transformed into a serial optimal flight mission planning problems through designing weighted combined cost functions, and the optimal mission planning problem that includes multi-thrust-coast phases can be transformed and solved rapidly by high precision pseudospectral knotting- nonlinear programming method. Numerical simulation shows that this method can generate controllable optimal flight trajectories while evaluating the integrated maneuver capability under various constraints including fuel consume and control angle rate limit, which will enhance the mission adaptability of finite thrust long-term in orbit maneuver vehicles.

Orbit maneuver; Capability envelop; Mission planning; Pseudospectral knotting

V412.4

A

1004-7182(2017)01-0014-04

10.7654/j.issn.1004-7182.20170104

2015-09-16；

2015-11-17

解永鋒（1984-），男，博士，工程師，主要研究方向?yàn)轱w行器動(dòng)力學(xué)控制