柯森锎，李爽，肖東東，王衛(wèi)華，聶欽博

（1. 南京航空航天大學(xué) 航天學(xué)院，南京 210016；2. 南京航空航天大學(xué) 航天新技術(shù)實(shí)驗(yàn)室，南京 210016；3. 上海航天控制技術(shù)研究所，上海 201109）

0 引 言

火星是離地球最近的類地行星，對(duì)火星的探測(cè)有助于人類進(jìn)一步理解認(rèn)識(shí)太陽(yáng)系的起源和演變，將宇宙科學(xué)推上新臺(tái)階[1]。在人類成功實(shí)施火星著陸探測(cè)之后，采樣返回將是火星探測(cè)的下一個(gè)里程碑[2-3]?；鹦巧仙鳎∕ars Ascent Vehicle，MAV）發(fā)射總質(zhì)量是關(guān)系到火星采樣返回任務(wù)可行性的關(guān)鍵因素，不僅關(guān)系到所能攜帶的火星樣本的質(zhì)量，而且直接影響火星表面上升軌跡設(shè)計(jì)。同時(shí)，它也關(guān)系到火星發(fā)射支持系統(tǒng)的質(zhì)量，進(jìn)一步影響到火星登陸器的質(zhì)量，最后傳導(dǎo)到地球上火箭的運(yùn)載能力要求上。MAV每增加1 kg，則地球發(fā)射質(zhì)量將相應(yīng)增加約8 kg[4]。而在當(dāng)前技術(shù)條件下每將1 kg的有效設(shè)備送上火星需要花費(fèi)約100萬(wàn)美元的成本[5]。所以MAV的設(shè)計(jì)不僅需要考慮火箭的運(yùn)載能力、火星進(jìn)入下降著陸（Entry，Descent and Landing，EDL）能力的可行性約束[6]，還面臨著經(jīng)濟(jì)性的考慮。

為了實(shí)現(xiàn)最輕的火星上升器設(shè)計(jì)，必須對(duì)MAV分級(jí)參數(shù)和上升軌跡進(jìn)行優(yōu)化。上升器參數(shù)包括各級(jí)火箭的干重、燃料質(zhì)量、推力、比沖等。軌跡優(yōu)化的經(jīng)典應(yīng)用領(lǐng)域?yàn)榛鸺淖顑?yōu)控制問(wèn)題[7-10]，在該問(wèn)題中給定火箭參數(shù)，要求以燃耗最優(yōu)為目標(biāo)函數(shù)進(jìn)行軌跡優(yōu)化以得到最優(yōu)的控制量。Betts對(duì)軌跡優(yōu)化問(wèn)題的求解方法進(jìn)行了概括分析，將其分為解析法和數(shù)值法，其中數(shù)值法又可以分為直接法和間接法[11]。高斯偽譜法是直接法的一種，Patterson提出了一種自適應(yīng)的高斯偽譜優(yōu)化方法并開(kāi)發(fā)了相應(yīng)的Matlab軟件GPOPS-II[12]。國(guó)內(nèi)的楊希祥基于高斯偽譜法做了固體火箭快速軌跡優(yōu)化研究[13]。如果火箭分級(jí)參數(shù)設(shè)計(jì)得不合理，那么軌跡優(yōu)化問(wèn)題就無(wú)法收斂得到理想的優(yōu)化結(jié)果，針對(duì)這一問(wèn)題Coskun提出了一種基于高斯偽譜法的優(yōu)化方法，能夠?qū)⒒鸺旨?jí)優(yōu)化和軌跡優(yōu)化耦合[14]。Benito也將高斯偽譜法引入到MAV的上升軌跡優(yōu)化問(wèn)題中來(lái)[15]。

本文以火星上升器總發(fā)射質(zhì)量最小為目標(biāo)函數(shù)，基于高斯偽譜法設(shè)計(jì)了一種針對(duì)兩級(jí)MAV的分級(jí)與軌跡耦合多階段優(yōu)化算法，最終得到MAV總發(fā)射質(zhì)量最小的火星上升器分級(jí)參數(shù)以及一條燃耗最優(yōu)上升軌跡。

1 問(wèn)題描述

1.1 飛行方案

在火星表面發(fā)射入軌所需的速度增量與地球相比大大減小，適合采用兩級(jí)MAV設(shè)計(jì)。本論文中上升器從火星表面起飛入軌的飛行方案如圖1所示，第1級(jí)點(diǎn)火開(kāi)始是一個(gè)很短的大推力的動(dòng)力上升段；然后是長(zhǎng)時(shí)間的無(wú)動(dòng)力滑行段一直到上升軌道的遠(yuǎn)拱點(diǎn)處；進(jìn)一步在遠(yuǎn)拱點(diǎn)處第2級(jí)火箭點(diǎn)火進(jìn)行軌道圓化進(jìn)入目標(biāo)圓軌道，然后釋放火星樣本艙，最終由火星軌道器實(shí)現(xiàn)對(duì)火星樣本艙的交匯捕獲以及返回地球[16-18]。

圖1 火星上升器飛行方案Fig. 1 MAV fly profile

1.2 坐標(biāo)系定義

1.3 兩級(jí)MAV質(zhì)量模型

一個(gè)一般化的兩級(jí)MAV的質(zhì)量模型如圖2所示。每一級(jí)的質(zhì)量由干重和燃料質(zhì)量組成[19]。而每級(jí)干重又由固定質(zhì)量和與燃料質(zhì)量相關(guān)的可變質(zhì)量組成。

圖2 一般化的兩級(jí)MAV質(zhì)量模型Fig. 2 The universal model for 2-stage MAVs

燃料質(zhì)量由每級(jí)的干重、載荷質(zhì)量、速度增量以及燃料比沖決定，燃料質(zhì)量為

第1級(jí)的載荷包括第2級(jí)總質(zhì)量加上有效載荷的質(zhì)量。

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

假設(shè)火星為旋轉(zhuǎn)的標(biāo)準(zhǔn)圓球，只考慮空氣阻力作用，在火心慣性直角坐標(biāo)系下建立MAV上升段的動(dòng)力學(xué)方程為

2 兩級(jí)MAV的分級(jí)與軌跡耦合多階段優(yōu)化算法

將兩級(jí)MAV的總發(fā)射質(zhì)量最小化問(wèn)題轉(zhuǎn)化為兩級(jí)MAV的分級(jí)與軌跡耦合多階段優(yōu)化問(wèn)題。根據(jù)兩級(jí)MAV的飛行方案將兩級(jí)MAV的分級(jí)與軌跡耦合優(yōu)化問(wèn)題分為3個(gè)階段：動(dòng)力段、滑行段、圓化段。

2.1 目標(biāo)函數(shù)與優(yōu)化變量

對(duì)于MAV的總發(fā)射質(zhì)量最小化問(wèn)題，其目標(biāo)是使得總發(fā)射質(zhì)量最小，等效于最小化第一階段的初始質(zhì)量，其目標(biāo)函數(shù)為

2.2 動(dòng)力學(xué)約束

動(dòng)力學(xué)約束由第1節(jié)中的動(dòng)力學(xué)模型給出。由于MAV總發(fā)射質(zhì)量最小化問(wèn)題需要對(duì)各級(jí)燃料質(zhì)量進(jìn)行優(yōu)化，所以增加了消耗的燃料質(zhì)量作為狀態(tài)量，需要增加式（10）所示的微分方程作為動(dòng)力學(xué)約束。

其中：p為消耗的燃料質(zhì)量。

2.3 邊界約束

本節(jié)將給出MAV總發(fā)射質(zhì)量最小化問(wèn)題各階段的邊界約束條件。

1）第1級(jí)推力段

2）無(wú)動(dòng)力滑行段

無(wú)動(dòng)力滑行段的初始條件由第1級(jí)推力段的終端狀態(tài)給出，即由2個(gè)階段的連接點(diǎn)約束決定。

為了增加MAV進(jìn)入火星穩(wěn)定軌道的概率，MAV將一直滑行到滑行軌道的遠(yuǎn)拱點(diǎn)再進(jìn)行軌道圓化，即該階段的終端約束為徑向速度為零，即

3）軌道圓化

該初始條件由滑行段的終端狀態(tài)給出，由兩個(gè)階段的連接點(diǎn)約束決定。在此階段的開(kāi)端MAV進(jìn)行一、二級(jí)火箭分離，拋棄第1級(jí)的同時(shí)第2級(jí)火箭點(diǎn)火進(jìn)行軌道圓化，所以第2級(jí)軌道圓化段與無(wú)動(dòng)力滑行段間的連接點(diǎn)約束要考慮第1級(jí)干重質(zhì)量。連接點(diǎn)約束為

將第2級(jí)軌道圓化段的初始燃料消耗質(zhì)量設(shè)置為零，即

該階段的終端條件為進(jìn)入目標(biāo)軌道，即滿足給定的軌道要素，與軌跡優(yōu)化問(wèn)題類似，但是需要增加終端質(zhì)量約束為

升交點(diǎn)赤經(jīng)可以通過(guò)改變發(fā)射時(shí)間來(lái)調(diào)整，因此本文中不考慮升交點(diǎn)赤經(jīng)約束，目標(biāo)軌道為圓軌道所以為目標(biāo)軌道高度。雖然圓軌道是數(shù)學(xué)上最簡(jiǎn)單的軌道，但是在現(xiàn)實(shí)中由于各種擾動(dòng)很難實(shí)現(xiàn)精確的圓軌道，如果將離心率設(shè)為零將會(huì)帶來(lái)一些數(shù)值計(jì)算上的困難，比如偽震蕩和不收斂。對(duì)于引入誤差范圍如下式所示。

2.4 路徑約束

1）第1級(jí)推力段

飛行高度約束

動(dòng)壓約束

推力角約束如式（31），推力角為推力方向與相對(duì)火星速度的夾角，其計(jì)算方程如式（32）所示，由于使用3DoF模型無(wú)法表示攻角，用推力角近似攻角。

2）無(wú)動(dòng)力滑行段

該階段的飛行高度約束以及動(dòng)壓約束與第1級(jí)推力段相同，因?yàn)檎麄€(gè)過(guò)程沒(méi)有推力，所以并沒(méi)有推力方向和推力角約束。

3）軌道圓化段

該階段在大氣層外所以無(wú)需考慮動(dòng)壓和推力角約束，但推力方向和飛行高度約束與第1級(jí)推力段相同。

3 高斯偽譜法

高斯偽譜法屬于直接法中的配點(diǎn)法，在處理分級(jí)與軌跡耦合優(yōu)化問(wèn)題時(shí)，通過(guò)在一系列拉格朗日高斯（Legendre Gauss，LG）點(diǎn)上構(gòu)建拉格朗日多項(xiàng)式來(lái)逼近狀態(tài)變量和控制變量[13]，進(jìn)一步對(duì)得到的多項(xiàng)式求導(dǎo)來(lái)近似狀態(tài)量的導(dǎo)數(shù)，將微分方程轉(zhuǎn)換為代數(shù)方程。對(duì)于終端條件和目標(biāo)函數(shù)中的積分項(xiàng)使用高斯積分來(lái)替換，將積分方程也轉(zhuǎn)換為代數(shù)方程，把分級(jí)與軌跡耦合優(yōu)化問(wèn)題離散后轉(zhuǎn)換為非線性規(guī)劃問(wèn)題[進(jìn)行]求解。

對(duì)（34）求導(dǎo)可以得到狀態(tài)量關(guān)于時(shí)間的導(dǎo)數(shù)為

基于式（36）可以將LG點(diǎn)處的動(dòng)力學(xué)方程約束轉(zhuǎn)換為代數(shù)方程形式為

終端狀態(tài)也需要滿足動(dòng)力學(xué)方程約束，將動(dòng)力學(xué)方程從初始時(shí)刻積分至終端時(shí)刻可以得到終端狀態(tài)量，使用高斯積分近似積分項(xiàng)，可以得到終端狀態(tài)量為

加上各離散點(diǎn)處的邊界條件和路徑約束就可以將兩級(jí)MAV分級(jí)與軌跡耦合優(yōu)化問(wèn)題離散化為非線性規(guī)劃問(wèn)題，求得滿足約束條件的狀態(tài)量和控制量使得目標(biāo)函數(shù)最小。

4 仿真實(shí)驗(yàn)與分析

對(duì)提出的基于高斯偽譜法的兩級(jí)MAV的分級(jí)與軌跡耦合多階段優(yōu)化算法仿真，目標(biāo)軌道參數(shù)、發(fā)射點(diǎn)位置以及兩級(jí)MAV質(zhì)量模型參數(shù)的詳細(xì)輸入值都由表1給出[19]。仿真程序在Matlab環(huán)境下運(yùn)行，計(jì)算機(jī)CPU為2.67 GHz Core i5，仿真計(jì)算得到兩級(jí)MAV的分級(jí)優(yōu)化結(jié)果如表2所示，同時(shí)得到一條標(biāo)稱上升軌跡如圖3～9所示。

表1 兩級(jí)MAV分級(jí)與軌跡耦合優(yōu)化算法仿真輸入?yún)?shù)Table 1 Input parameters of two-stage MAV coupled stagingtrajectory optimization algorithm simulation

表2 兩級(jí)MAV分級(jí)優(yōu)化結(jié)果Table 2 The result of two-stage MAV staging optimization

圖3 上升高度隨時(shí)間變化曲線Fig. 3 Altitude vs. time

圖4 飛行速度隨時(shí)間變化曲線Fig. 4 Speed vs. time

圖5 軌道傾角隨時(shí)間變化曲線Fig. 5 Inclination vs. time

由圖3～5可知滑行段時(shí)間約為565 s，燃耗最優(yōu)的軌跡并不是在第1級(jí)耗盡時(shí)達(dá)到目標(biāo)軌道傾角，而是留下約2.5°的軌道傾角誤差留給第2級(jí)進(jìn)行調(diào)整；由圖6～7第1級(jí)控制量可得發(fā)射方位角為42.5°，仰角為40.5°，第2級(jí)進(jìn)行軌道圓化的推力方向接近于固定方向；由于MAV的氣動(dòng)阻力系數(shù)Cd在這里被設(shè)定為簡(jiǎn)單的恒定值0.2[21]，所以得到的動(dòng)壓值只能作為一個(gè)概略的定性值，由圖8知由于燃燒時(shí)間短，第1級(jí)耗盡點(diǎn)高度較低，在第1級(jí)耗盡點(diǎn)處達(dá)到了動(dòng)壓的峰值，可以考慮增加燃燒時(shí)間以提高耗盡點(diǎn)高度，以降低耗盡點(diǎn)的動(dòng)壓使其為在峰值之后的一個(gè)較小值，可以減小滑行段姿態(tài)保持的難度；圖9中第1級(jí)的推力角起始值約為60°，之后快速減小并保持在5°以下，這是因?yàn)橥屏鞘峭屏εcMAV相對(duì)火星速度的夾角，而在起始時(shí)刻MAV相對(duì)火星速度為零，這里的起始推力角是由數(shù)值計(jì)算誤差得到的。

圖6 第一階段的控制量隨時(shí)間變化曲線Fig. 6 Control vs. time for phase 1

圖7 第三階段的控制量隨時(shí)間變化曲線Fig. 7 Control vs. time for phase 2

圖8 動(dòng)壓隨時(shí)間變化曲線Fig. 8 Dynamic pressure vs. time

圖9 第1級(jí)推力角隨時(shí)間變化曲線Fig. 9 Thrust angle vs. time for phase 1

通過(guò)仿真計(jì)算可知，該優(yōu)化算法的初值在給定范圍內(nèi)隨機(jī)選取都能在約7 min內(nèi)收斂得到最優(yōu)結(jié)果，計(jì)算精度達(dá)到10-6量級(jí)。

5 結(jié) 論

以火星采樣返回任務(wù)中的兩級(jí)火星上升器為研究對(duì)象，利用高斯偽譜法設(shè)計(jì)了在復(fù)雜多約束條件下兩級(jí)MAV分級(jí)與軌跡耦合多階段優(yōu)化算法，以求得總發(fā)射質(zhì)量最小的兩級(jí)MAV分級(jí)設(shè)計(jì)以及燃耗最優(yōu)的上升軌跡。該算法將MAV的分級(jí)優(yōu)化與軌跡優(yōu)化耦合在一起，解決了由于不合理的分級(jí)參數(shù)導(dǎo)致軌跡優(yōu)化無(wú)法收斂的問(wèn)題。通過(guò)仿真驗(yàn)證可以得到如下結(jié)論。

1）為達(dá)到燃耗最優(yōu)的目的，兩級(jí)MAV的第1級(jí)將留下一定的軌道傾角偏轉(zhuǎn)量由第2級(jí)完成。

2）為減小兩級(jí)MAV在滑行段的姿控難度可以適當(dāng)減小第1級(jí)火箭的推力進(jìn)而減小第1級(jí)耗盡點(diǎn)的動(dòng)壓。

3）提出的優(yōu)化算法對(duì)于初值選擇的敏感度較小，具有較強(qiáng)的魯棒性，且收斂速度較快，方便于進(jìn)行快速的兩級(jí)MAV設(shè)計(jì)優(yōu)化。

本文僅考慮了總發(fā)射質(zhì)量最小和燃耗最優(yōu)作為性能指標(biāo)的火星表面上升段軌跡設(shè)計(jì)問(wèn)題，而在實(shí)際應(yīng)用中還需要考慮飛行時(shí)間、火星軌道器觀測(cè)跟蹤范圍等其他因素，未來(lái)需要針對(duì)MAV上升軌跡的多性能指標(biāo)優(yōu)化開(kāi)展進(jìn)一步的研究。