魯 鵬，李雅軒，劉新福

（北京理工大學(xué)宇航學(xué)院，北京，100081）

0 引 言

航天運(yùn)輸系統(tǒng)是開(kāi)展空間活動(dòng)的核心支撐，是發(fā)展航天事業(yè)、建設(shè)航天強(qiáng)國(guó)的基本前提，是加速發(fā)展空間科學(xué)技術(shù)、和平開(kāi)發(fā)空間資源和維護(hù)國(guó)家空間安全的重要保障［1］。近年來(lái)，隨著人類(lèi)空間資源開(kāi)發(fā)需求的快速增長(zhǎng)，可重復(fù)使用運(yùn)載器受到了全球的廣泛關(guān)注，世界各主要航天大國(guó)積極開(kāi)展了相關(guān)的理論研究和工程實(shí)踐［2］。

瞄準(zhǔn)星際移民和深空探測(cè)的目標(biāo)，美國(guó)太空探索技術(shù)公司（SpaceX）研制了名為星艦（Starship）的新一代可重復(fù)使用航天運(yùn)輸系統(tǒng)。2016 年9 月，SpaceX 公司正式對(duì)外公布了“星際運(yùn)輸系統(tǒng)”方案［3］。此后，星艦設(shè)計(jì)方案經(jīng)歷了多次重大改進(jìn)和演化。自2020 年起，通過(guò)開(kāi)展密集的飛行測(cè)試，SpaceX 公司逐步掌握了150 m 低空懸停、10 km 高空飛行、翻轉(zhuǎn)機(jī)動(dòng)和垂直定點(diǎn)軟著陸等關(guān)鍵技術(shù)，項(xiàng)目取得重大進(jìn)展。值得關(guān)注的是，星艦采取了與傳統(tǒng)航天器截然不同的著陸方式：著陸前，星艦采用90°左右大攻角飛行，在視覺(jué)上以腹部平拍的姿態(tài)向地面飛行；接近地面時(shí)，星艦發(fā)動(dòng)機(jī)開(kāi)機(jī)進(jìn)行推力矢量控制，聯(lián)合首尾四個(gè)氣動(dòng)舵，完成姿態(tài)由水平向垂直的快速翻轉(zhuǎn)，并控制星艦減速和修正位置誤差，實(shí)現(xiàn)定點(diǎn)軟著陸［4］。圖1展示了該特殊的著陸過(guò)程。這種獨(dú)特的著陸方式對(duì)飛行器的軌跡優(yōu)化設(shè)計(jì)提出了新的挑戰(zhàn)。大范圍的姿態(tài)機(jī)動(dòng)和復(fù)雜的約束條件必須在優(yōu)化設(shè)計(jì)翻轉(zhuǎn)著陸飛行軌跡時(shí)得到充分的考慮。

圖1 類(lèi)星艦飛行器翻轉(zhuǎn)著陸過(guò)程示意Fig.1 The flipping and landing process of a starship-like vehicle

近年來(lái)，凸優(yōu)化方法被廣泛地應(yīng)用于飛行器軌跡優(yōu)化設(shè)計(jì)，特別是在解決星際探測(cè)器、運(yùn)載火箭等常規(guī)航天器的著陸軌跡優(yōu)化問(wèn)題時(shí)發(fā)揮積極作用［5］。美國(guó)噴氣推進(jìn)實(shí)驗(yàn)室的A??kme?e 等［6］通過(guò)使用無(wú)損松弛和變量替換等方法將帶有非凸推力約束的火星著陸軌跡優(yōu)化問(wèn)題轉(zhuǎn)化為凸優(yōu)化問(wèn)題，顯著提高了軌跡優(yōu)化算法的可靠性和計(jì)算效率。Szmuk等［7］對(duì)六自由度火星動(dòng)力下降制導(dǎo)問(wèn)題進(jìn)行了深入研究，通過(guò)設(shè)計(jì)序列凸優(yōu)化算法計(jì)算考慮姿態(tài)動(dòng)力學(xué)的著陸飛行軌跡。相比于火星著陸，因?yàn)闅鈩?dòng)力不可忽略，運(yùn)載火箭著陸在地球表面的軌跡優(yōu)化問(wèn)題求解難度更大。劉新福［8］將推力和攻角同時(shí)作為控制量，把對(duì)應(yīng)的燃料最優(yōu)著陸問(wèn)題順利轉(zhuǎn)化為凸優(yōu)化問(wèn)題進(jìn)行迭代求解，并理論證明了所采用松弛技術(shù)的有效性。楊潤(rùn)秋等［9］將該研究工作拓展到飛行時(shí)間自由的問(wèn)題，用高度代替時(shí)間作為自變量對(duì)優(yōu)化問(wèn)題進(jìn)行重構(gòu)，并基于序列凸優(yōu)化設(shè)計(jì)了軌跡優(yōu)化算法，其中使用的非線(xiàn)性保留技術(shù)被證實(shí)可以有效提升算法收斂性。

盡管上述方法能夠有效解決常規(guī)航天器的著陸軌跡優(yōu)化問(wèn)題，但對(duì)類(lèi)星艦飛行器的翻轉(zhuǎn)著陸軌跡優(yōu)化并不太適用。原因在于類(lèi)星艦飛行器在滑翔段和著陸段之間存在獨(dú)特的翻轉(zhuǎn)過(guò)程。大范圍的姿態(tài)機(jī)動(dòng)導(dǎo)致類(lèi)星艦飛行器不能被簡(jiǎn)單地處理為質(zhì)點(diǎn)模型，而必須在設(shè)計(jì)著陸軌跡時(shí)考慮姿態(tài)動(dòng)力學(xué)和因姿態(tài)控制產(chǎn)生的復(fù)雜約束條件。但是，目前針對(duì)此類(lèi)軌跡優(yōu)化問(wèn)題的研究非常有限。

本文針對(duì)類(lèi)星艦飛行器的翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題進(jìn)行了深入研究。首先根據(jù)類(lèi)星艦飛行器在翻轉(zhuǎn)著陸階段的飛行力學(xué)特征和任務(wù)特點(diǎn)，建立考慮姿態(tài)機(jī)動(dòng)和典型約束條件的翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題。然后研究將該優(yōu)化問(wèn)題通過(guò)凸化方法轉(zhuǎn)化為凸優(yōu)化問(wèn)題，進(jìn)而設(shè)計(jì)迭代算法序列求解。通過(guò)提供合理的初值、引入虛擬控制和在目標(biāo)函數(shù)中懲罰信賴(lài)域半徑等方法，該迭代算法表現(xiàn)出良好的收斂性。針對(duì)發(fā)動(dòng)機(jī)擺角振蕩問(wèn)題，本文提出將發(fā)動(dòng)機(jī)擺角響應(yīng)引入動(dòng)力學(xué)，這種處理被證實(shí)能夠有效抑制發(fā)動(dòng)機(jī)擺角振蕩，并有利于迭代算法的收斂。此外，依托軌跡優(yōu)化算法，本文亦分析了發(fā)動(dòng)機(jī)個(gè)數(shù)和優(yōu)化目標(biāo)選取等因素對(duì)類(lèi)星艦飛行器翻轉(zhuǎn)著陸過(guò)程的影響，相關(guān)數(shù)值結(jié)果可以為類(lèi)星艦飛行器的翻轉(zhuǎn)著陸軌跡優(yōu)化設(shè)計(jì)提供參考。

1 翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題

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

本文將研究類(lèi)星艦飛行器在二維平面內(nèi)考慮一個(gè)姿態(tài)角的翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題。再入過(guò)程中飛行器控制能力強(qiáng)，通?？墒癸w行器再入終端速度矢量和著陸點(diǎn)位于同一垂直平面內(nèi)，而著陸過(guò)程中的側(cè)向運(yùn)動(dòng)通常會(huì)消耗額外燃料，因此在進(jìn)行類(lèi)星艦飛行器軌跡優(yōu)化時(shí)可以忽略側(cè)向運(yùn)動(dòng)。實(shí)際飛行時(shí)，著陸過(guò)程中的側(cè)向位置誤差可以通過(guò)控制進(jìn)行修正。構(gòu)造類(lèi)星艦飛行器翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題時(shí)會(huì)使用到著陸坐標(biāo)系Oxy和體坐標(biāo)系Obxbyb，如圖2所示。

圖2 著陸坐標(biāo)系Oxy和體坐標(biāo)系Obxb ybFig.2 The definition of landing coordinate system Oxy and the body-fixed coordinate system Obxb yb.

圖3a 展示了飛行器翻轉(zhuǎn)著陸過(guò)程受力分析結(jié)果，根據(jù)牛頓第二定律建立動(dòng)力學(xué)方程。為了避免在問(wèn)題求解時(shí)出現(xiàn)數(shù)值問(wèn)題，需要對(duì)動(dòng)力學(xué)方程中的物理量進(jìn)行無(wú)量綱化處理。長(zhǎng)度的無(wú)量綱系數(shù)為飛行器初始高度ry0，速度無(wú)量綱系數(shù)為ry0g0（g0是地球表面重力加速度），時(shí)間無(wú)量綱系數(shù)為ry0g0，質(zhì)量無(wú)量綱系數(shù)為初始質(zhì)量m0。著陸坐標(biāo)系下，無(wú)量綱的飛行器翻轉(zhuǎn)著陸的動(dòng)力學(xué)方程，記為x?(t) =f(x(t)，u(t))，如下：

圖3 飛行器受力分析和結(jié)構(gòu)參數(shù)Fig.3 Force analysis and structural parameters of the vehicle

式中r∈R3為飛行器位置矢量；v∈R3為飛行器速度矢量；m為飛行器質(zhì)量；T為發(fā)動(dòng)機(jī)推力矢量；T為推力矢量的模，T=‖T‖；D為氣動(dòng)阻力矢量；g為重力加速度矢量；θ為飛行器俯仰角（定義見(jiàn)圖2）；ω為俯仰角速度；MT為推力產(chǎn)生的力矩；MD為阻力力矩；Jz為飛行器繞質(zhì)心的轉(zhuǎn)動(dòng)慣量；αe為質(zhì)量消耗率，αe=-1Isp，Isp為發(fā)動(dòng)機(jī)比沖；δ為發(fā)動(dòng)機(jī)擺角；δd為δ的響應(yīng)；Td為一階慣性環(huán)節(jié)的無(wú)量綱時(shí)間常數(shù)。本文將發(fā)動(dòng)機(jī)擺角視為控制變量，它同時(shí)控制位置和姿態(tài)，進(jìn)行軌跡優(yōu)化時(shí)會(huì)導(dǎo)致嚴(yán)重的發(fā)動(dòng)機(jī)擺角振蕩。鑒于此，本文將發(fā)動(dòng)機(jī)擺角響應(yīng)引入動(dòng)力學(xué)（見(jiàn)式（6））。數(shù)值仿真結(jié)果表明，這一思路將有效抑制發(fā)動(dòng)機(jī)擺角振蕩現(xiàn)象，進(jìn)而提高了2.3 節(jié)設(shè)計(jì)的軌跡優(yōu)化算法的收斂性。x為狀態(tài)變量，x=[rTvTθ ω m δd]T，u為控制變量，u=[T δ]T。本文將動(dòng)力學(xué)投影在著陸坐標(biāo)系中，因此無(wú)量綱的推力T、阻力D以及相應(yīng)推力力矩MT和阻力力矩MD的計(jì)算公式如下：

1.2 其它約束條件

飛行器翻轉(zhuǎn)著陸軌跡優(yōu)化不僅受到動(dòng)力學(xué)方程的約束，還有端點(diǎn)約束和過(guò)程約束，接下來(lái)詳細(xì)介紹這些約束。

飛行器初始狀態(tài)約束可以表述成Φ= 0 的形式。本優(yōu)化問(wèn)題給定初始位置r(1)0、速度v(1)0、俯仰角θ(1)0、俯仰角速度ω(1)0和初始質(zhì)量m(1)0，本文使用上標(biāo)“（1）”標(biāo)記翻轉(zhuǎn)段對(duì)應(yīng)的量，使用上標(biāo)“（2）”標(biāo)記著陸段對(duì)應(yīng)的量。因此Φ具體表達(dá)式如下：

翻轉(zhuǎn)著陸過(guò)程中飛行器不能有向上的速度，即豎直速度vy不能為正，因此得到速度約束：

為了保證翻轉(zhuǎn)過(guò)程中飛行器姿態(tài)在合理范圍，保證著陸過(guò)程飛行器垂直著陸，翻轉(zhuǎn)著陸過(guò)程施加如下關(guān)于俯仰角的邊界約束：

考慮到飛行器在載人時(shí)的舒適性和姿態(tài)控制系統(tǒng)的能力，姿態(tài)角速度不能超過(guò)最大值ωmax，即有：

發(fā)動(dòng)機(jī)的最大擺角記為δmax，故擺角需要滿(mǎn)足：

發(fā)動(dòng)機(jī)擺角速度不能超過(guò)最大值δ?max，可以通過(guò)約束發(fā)動(dòng)機(jī)擺角響應(yīng)的變化率δ?d實(shí)現(xiàn)對(duì)發(fā)動(dòng)機(jī)擺角速度的限制，由式（6）可知，需引入如下限制發(fā)動(dòng)機(jī)擺角速度的約束：

式中n為發(fā)動(dòng)機(jī)個(gè)數(shù)，n∈N+。

為了方便垂直著陸，每個(gè)時(shí)刻飛行器和著陸點(diǎn)連線(xiàn)與豎直方向的夾角不能大于給定角度γgs，這就需要添加下滑道約束：

式中ry為位置矢量在著陸坐標(biāo)系y軸方向的分量。該下滑道約束是一個(gè)二階錐約束。

過(guò)程約束（12）～（18）可以改寫(xiě)為如下不等式約束：

末端等式約束可以表示為E=0 的形式，本優(yōu)化問(wèn)題中飛行器將著陸到給定位置r(2)f，與此同時(shí)末端速度到達(dá)給定值v(2)f、俯仰角到達(dá)θ(2)f并且俯仰角速度到達(dá)ω(2)f。因此，E具有如下形式：

由于燃料有限，終端質(zhì)量需要滿(mǎn)足m(2)(t(2)f)≥mdry，mdry為燃料耗盡時(shí)飛行器的質(zhì)量。

飛行器在翻轉(zhuǎn)段和著陸段交接點(diǎn)處的飛行狀態(tài)應(yīng)當(dāng)連續(xù)，連接約束可以表示為Ψ=0，其中：

翻轉(zhuǎn)段耗時(shí)、著陸段耗時(shí)和翻轉(zhuǎn)著陸段總飛行時(shí)間滿(mǎn)足約束：

式中tmin為總飛行時(shí)間下界；tmax為總飛行時(shí)間上界。

1.3 翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題

在優(yōu)化類(lèi)星艦飛行器翻轉(zhuǎn)著陸的軌跡時(shí)，為節(jié)約燃料，同時(shí)保證翻轉(zhuǎn)結(jié)束時(shí)飛行器高度不能太低，設(shè)計(jì)了如下目標(biāo)函數(shù)：

式中μ為懲罰系數(shù)，調(diào)節(jié)μ可以改變翻轉(zhuǎn)段結(jié)束時(shí)飛行器的離地高度。當(dāng)μ= 0 時(shí)，該優(yōu)化目標(biāo)將轉(zhuǎn)化為翻轉(zhuǎn)著陸過(guò)程燃料最優(yōu)。本文3.4 節(jié)將詳細(xì)探討該懲罰系數(shù)的取值對(duì)翻轉(zhuǎn)著陸軌跡的影響。

綜上所述，可以構(gòu)造翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題P0：

問(wèn)題P0的本質(zhì)是多階段飛行時(shí)間自由的最優(yōu)控制問(wèn)題。因?yàn)閯?dòng)力學(xué)約束非線(xiàn)性很強(qiáng)，P0 求解難度較大。

2 凸優(yōu)化求解

本節(jié)將介紹使用改進(jìn)的序列凸優(yōu)化（Sequential Convex Programming，SCP）算法求解飛行時(shí)間自由的類(lèi)星艦飛行器翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題的方法。雖然序列凸優(yōu)化算法無(wú)法保證收斂，但是可以使用一些方法提高其收斂性。為了提高收斂性，求解問(wèn)題P0時(shí)引入虛擬控制，懲罰信賴(lài)域半徑［10］，并給出了初值選取方法。為了公式的簡(jiǎn)潔性，本節(jié)不再添加區(qū)分翻轉(zhuǎn)段和著陸段的上標(biāo)“（i）”。

為了將飛行時(shí)間自由的最優(yōu)控制問(wèn)題P0轉(zhuǎn)化為等價(jià)的飛行時(shí)間固定的最優(yōu)控制問(wèn)題，引入變量τ=(t-t0)/Δt，其中Δt=tf-t0表示該階段耗時(shí)，t0是該階段的初始時(shí)間，tf是該階段的終端時(shí)間。將動(dòng)力學(xué)方程轉(zhuǎn)換為以τ∈[0，1]為自變量的形式：

式中 上標(biāo)“'”表示變量對(duì)τ求導(dǎo)。狀態(tài)變量和控制變量是t的函數(shù)，同樣也是τ的函數(shù)，故有：

當(dāng)以τ作為自變量時(shí)，初始狀態(tài)約束是τ= 0時(shí)的狀態(tài)約束，終端狀態(tài)約束是τ= 1時(shí)的狀態(tài)約束。

2.1 線(xiàn)性化

線(xiàn)性化可以將上述非線(xiàn)性動(dòng)力學(xué)方程（25）轉(zhuǎn)化為線(xiàn)性動(dòng)力學(xué)方程。將方程（25）在上一次迭代得到的解附近泰勒展開(kāi)并保留常數(shù)項(xiàng)和一階項(xiàng)，用展開(kāi)后的線(xiàn)性方程近似原本的非線(xiàn)性動(dòng)力學(xué)方程。這里將上一次迭代得到的解記為z[k](τ) ?[(x[k])T(u[k])T]T，時(shí)間記為Δt[k]，那么在z[k](τ)和Δt[k]附近線(xiàn)性化方程（25），可以得到如下線(xiàn)性動(dòng)力學(xué)方程：

為了構(gòu)造凸優(yōu)化問(wèn)題，除了動(dòng)力學(xué)方程外其它非凸約束也需要進(jìn)行線(xiàn)性化處理，例如：假設(shè)sn(z)是關(guān)于z的非凸函數(shù)，則如下約束是非凸約束：

使用一階泰勒展開(kāi)，略去高階項(xiàng)，sn(z(τ))≤0 可以由如下線(xiàn)性不等式約束近似：

而問(wèn)題P0 中除了動(dòng)力學(xué)約束外沒(méi)有其它非凸約束。

2.2 離 散

為了將連續(xù)時(shí)間最優(yōu)控制問(wèn)題轉(zhuǎn)化為非線(xiàn)性規(guī)劃問(wèn)題，需要對(duì)連續(xù)時(shí)間最優(yōu)控制問(wèn)題進(jìn)行離散。將τ所在的區(qū)間［0，1］離散成N份，0 =τ0<τ1< …<τN=1，式中τ0，τ1，…，τN是離散節(jié)點(diǎn)，與自變量對(duì)應(yīng)的狀態(tài) 離 散 點(diǎn) 記 為x0，x1，…，xN，控 制 離 散 點(diǎn) 記 為u0，u1，…，uN。離散動(dòng)力學(xué)的數(shù)值方法有很多［11］，包括歐拉法、梯形法、零階保持器和一階保持器等，本文使用梯形法。線(xiàn)性化后的動(dòng)力學(xué)方程（26）經(jīng)梯形法離散得到線(xiàn)性等式：

式中hj=τj+1-τj。

動(dòng)力學(xué)以外的其它約束，例如sn(z(τ)) ≤0，離散形式是每個(gè)離散節(jié)點(diǎn)上都滿(mǎn)足相應(yīng)約束：

2.3 虛擬控制和懲罰信賴(lài)域

非線(xiàn)性動(dòng)力學(xué)和其它非線(xiàn)性約束經(jīng)過(guò)線(xiàn)性化后，可能使原本可行的問(wèn)題變成不可行的問(wèn)題。為了解決由線(xiàn)性化導(dǎo)致的問(wèn)題不可行，在動(dòng)力學(xué)方程（29）上添加虛擬控制。離散的動(dòng)力學(xué)方程（29）添加虛擬控制vc(τ) ∈R8后 得 到 等 式 約 束（31）， 記 為

式中vcj=vc(τj)，虛擬控制vcj可以保證任意下一個(gè)節(jié)點(diǎn)處的狀態(tài)xj+1都是可以達(dá)到的。需要注意，如果序列凸優(yōu)化算法收斂時(shí)虛擬控制的值不為零，則得到的結(jié)果就不滿(mǎn)足線(xiàn)性化前的動(dòng)力學(xué)約束。為了使程序收斂時(shí)虛擬控制為零，需要在目標(biāo)函數(shù)中增加關(guān)于虛擬控制的懲罰項(xiàng)Jvc，并且需要為虛擬控制設(shè)置一個(gè)較大的懲罰系數(shù)λvc，本文使用的Jvc如下：

為了保證線(xiàn)性化的合理性，需要添加信賴(lài)域約束。本文使用可變信賴(lài)域，信賴(lài)域作為懲罰項(xiàng)放在目標(biāo)函數(shù)中，即在目標(biāo)函數(shù)中添加關(guān)于信賴(lài)域的懲罰項(xiàng)Jtr：

式中λtr為懲罰系數(shù)。這種以添加懲罰項(xiàng)的方式處理信賴(lài)域半徑的方法通?？梢允箚?wèn)題收斂得更快。λtr取值一般較小，因?yàn)楫?dāng)該值較大時(shí)，前后兩次迭代的結(jié)果相差較小，會(huì)導(dǎo)致序列凸優(yōu)算法收斂速度變慢，并且很容易收斂到初值附近的可行解。

至此，可以得到時(shí)間自由的兩階段軌跡優(yōu)化問(wèn)題P0被凸化后的優(yōu)化問(wèn)題P1。本文在翻轉(zhuǎn)段將τ離散為N1等份，在著陸段將τ離散為N2等份。

序列求解凸優(yōu)化問(wèn)題P1，當(dāng)Jvc< ?vc并且Jtr

圖4 軌跡優(yōu)化算法流程Fig.4 Flow chart of trajectory optimization algorithm

2.4 初值選取

對(duì)于動(dòng)力學(xué)和其它約束較為簡(jiǎn)單的優(yōu)化問(wèn)題，在使用序列凸優(yōu)化算法求解時(shí)，通?？梢灾苯邮褂镁€(xiàn)性初值。而對(duì)于本文構(gòu)造的動(dòng)力學(xué)非線(xiàn)性很強(qiáng)的優(yōu)化問(wèn)題，由于線(xiàn)性初值表現(xiàn)不佳，很容易造成算法1收斂慢或無(wú)法收斂的情況。為了提高算法1的收斂性，在這里給出一種選取更好初值的方法。

在求翻轉(zhuǎn)段初值時(shí)，先忽略翻轉(zhuǎn)過(guò)程中氣動(dòng)阻力產(chǎn)生的力矩，然后給定發(fā)動(dòng)機(jī)推力和擺角，積分到俯仰角等于給定的θs且垂直速度等于vys時(shí)結(jié)束，積分結(jié)果作為翻轉(zhuǎn)段初值。θs取90°附近的值便于后續(xù)垂直著陸。速度vys不能太大，需要保證著陸時(shí)速度能減到0。在計(jì)算著陸段初值時(shí)，求解在問(wèn)題P0著陸段基礎(chǔ)上簡(jiǎn)化得到的凸優(yōu)化問(wèn)題，以簡(jiǎn)化問(wèn)題的解作為著陸段初值。具體來(lái)說(shuō)就是給定簡(jiǎn)化問(wèn)題的著陸時(shí)間，忽略氣動(dòng)力，關(guān)于速度的動(dòng)力學(xué)方程不考慮質(zhì)量變化，不考慮姿態(tài)動(dòng)力學(xué)和發(fā)動(dòng)機(jī)擺角，并滿(mǎn)足部分過(guò)程和終端約束，即求解如下凸優(yōu)化問(wèn)題：

式中x0取翻轉(zhuǎn)段末端狀態(tài)；γT為推力方向和y軸間允許的最大夾角；ε為松弛變量；不使用ε計(jì)算初值的優(yōu)化問(wèn)題可能無(wú)解。 上述優(yōu)化問(wèn)題中不包括θ，ω，δ，δd。俯仰角θ根據(jù)推力方向計(jì)算，引入推力方向約束就是為了避免通過(guò)推力方向計(jì)算的俯仰角過(guò)大，ω通過(guò)對(duì)θ進(jìn)行差商得到，δ(t) =δd(t) =0。求解上述凸優(yōu)化問(wèn)題可以得到一個(gè)較好的初值，且計(jì)算效率高。仿真部分將比較本初值選取方法和線(xiàn)性初值選取方法對(duì)算法1的影響。

3 數(shù)值仿真

本節(jié)將通過(guò)數(shù)值仿真來(lái)分析類(lèi)星艦飛行器翻轉(zhuǎn)著陸問(wèn)題的特性，希望能為相關(guān)任務(wù)設(shè)計(jì)提供一定的參考。本節(jié)數(shù)值仿真主要參考星艦公開(kāi)的參數(shù)［13］，部分重要參數(shù)見(jiàn)表1。本文在實(shí)現(xiàn)算法1 時(shí)使用了凸優(yōu)化求解器ECOS［14］。

表1 飛行器參數(shù)Tab.1 Vehicle configurations

飛行器高l= 50 m，lcg= 0.6l，lcp= 0.55l。飛行器初始位置r0=[360 1200]Tm，初始速度v0=［-18.82-106.73]Tm/s，θ0= 170°，ω0= 0（°）/s，m0= 132 t。用于計(jì)算的初值θs= 80°，vys=-20 m/s。翻轉(zhuǎn)過(guò)程使用3 臺(tái)發(fā)動(dòng)機(jī)，著陸使用2 臺(tái)發(fā)動(dòng)機(jī)。末端位置rf=[0 0]Tm，末端速度vf=[0 -0.1]Tm/s，θf(wàn)= 90°，ωf=0（°）/s。邊界約束中θ(1)max= 170°，θ(1)min= 45°，θ(2)max=105°，θ(2)min= 75°。算 法1 參 數(shù)：N1= 40，N2= 50，λvc= 103，λtr= 5 × 10-5，收斂判據(jù)?vc= 1 × 10-8，?tr=1 × 10-8，目標(biāo)函數(shù)中μ=0.3。

使用圖4算法求解該算例時(shí)，迭代7次收斂，翻轉(zhuǎn)段耗時(shí)2.43 s，著陸段耗時(shí)13.61 s，整個(gè)過(guò)程消耗燃料9.773 t。圖5 展示了優(yōu)化結(jié)果：圖5a 展示了翻轉(zhuǎn)著陸過(guò)程的路徑點(diǎn)和推力方向。圖5b展示了俯仰角和俯仰角速度的變化，最大的俯仰角速度小于給定上界36 。優(yōu)化的控制量繪制在圖5c和圖5d中，由于目標(biāo)函數(shù)考慮了讓翻轉(zhuǎn)過(guò)程高度差較小，所以翻轉(zhuǎn)時(shí)先以最大推力和最大擺角加速翻轉(zhuǎn)，接著使用小推力翻轉(zhuǎn)，同時(shí)擺角開(kāi)始減小為接下來(lái)減小俯仰角速度做準(zhǔn)備，然后在著陸段開(kāi)始時(shí)以最大推力和最小擺角為翻轉(zhuǎn)過(guò)程減速，即降低俯仰角速度。整個(gè)翻轉(zhuǎn)著陸過(guò)程推力幾乎是砰—砰（bang-bang）形式，發(fā)動(dòng)機(jī)擺角不大于10°。

圖5 飛行器翻轉(zhuǎn)著陸軌跡優(yōu)化結(jié)果Fig.5 Optimized results of the flipping and landing trajectory of the Starship-like vehicle

接下來(lái)改變部分仿真條件，深入研究類(lèi)星艦飛行器翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題的性質(zhì)，并將對(duì)應(yīng)結(jié)果與圖5中的結(jié)果進(jìn)行比較。

3.1 不同初值選取方法的影響

對(duì)于復(fù)雜的軌跡優(yōu)化問(wèn)題，初值選取對(duì)于序列凸優(yōu)化算法的收斂性有較大影響。這一小節(jié)將比較線(xiàn)性初值和本文初值選取方法對(duì)算法1的影響。

使用線(xiàn)性初值，算法1 迭代9 次后收斂，比使用本文初值選取方法多了2次迭代，從圖6a可以看出使用本文初值收斂更快。翻轉(zhuǎn)段耗時(shí)2.38 s，著陸段耗時(shí)14.02 s，整個(gè)過(guò)程消耗燃料9.947 t，比使用本文初值多消耗174 kg。圖6c 展示了推力剖面的對(duì)比，使用線(xiàn)性初值時(shí)，推力不是砰—砰形式。此外，使用線(xiàn)性初值時(shí)，發(fā)動(dòng)機(jī)擺角剖面也不如使用本文初值時(shí)那么平緩，如圖6d所示。大量仿真實(shí)例表明算法1使用本文初值選取方法比線(xiàn)性初值收斂性更好，在有些使用線(xiàn)性初值無(wú)法收斂的情況下，使用本文初值也能收斂。本文初值選取方法不僅提高了算法的收斂性，而且優(yōu)化的推力和發(fā)動(dòng)機(jī)擺角剖面更適合工程使用。

圖6 使用線(xiàn)性初值和本文方法確定的初值的結(jié)果對(duì)比Fig.6 The comparison between the results with linear initial guess or the initial guess obtained by the proposed method

3.2 發(fā)動(dòng)機(jī)擺角響應(yīng)相關(guān)約束的作用

在3.1 節(jié)中說(shuō)明了引入發(fā)動(dòng)機(jī)擺角響應(yīng)和相關(guān)約束的作用——避免求解的發(fā)動(dòng)機(jī)擺角振蕩。這一小節(jié)將給出不使用發(fā)動(dòng)機(jī)擺角響應(yīng)δd而其它部分都相同時(shí)的計(jì)算結(jié)果。

如前所述，使用δd的算例迭代7 次收斂，而不使用δd的算例迭代20次未收斂。圖7給出了使用和不使用δd計(jì)算結(jié)果的對(duì)比（不使用擺角響應(yīng)的算例20次迭代未收斂，這里畫(huà)了第20次迭代中的解），圖7a和圖7c顯示兩個(gè)算例計(jì)算的推力和速度剖面幾乎相同，但是發(fā)動(dòng)機(jī)擺角差別較大，不使用δd的算例得到的發(fā)動(dòng)機(jī)擺角（見(jiàn)圖7b 中紅色實(shí)線(xiàn)）著陸前出現(xiàn)了劇烈的振蕩。從圖7d可以看出不使用δd時(shí)，Jtr還未小于?tr就不再減小，因此不使用δd的算例沒(méi)有收斂。發(fā)動(dòng)機(jī)擺角劇烈振蕩是此算例不好收斂的主要原因。

圖7 不使用和使用擺角響應(yīng)約束求解結(jié)果的對(duì)比Fig.7 The results solved with and without δd

3.3 著陸使用發(fā)動(dòng)機(jī)臺(tái)數(shù)的影響

仿真使用的飛行器的發(fā)動(dòng)機(jī)推力較大，足以實(shí)現(xiàn)單臺(tái)發(fā)動(dòng)機(jī)著陸。這一小節(jié)給出單臺(tái)發(fā)動(dòng)機(jī)和兩臺(tái)發(fā)動(dòng)機(jī)著陸時(shí)的對(duì)比結(jié)果。翻轉(zhuǎn)段仍然使用三臺(tái)發(fā)動(dòng)機(jī)。

使用單臺(tái)和兩臺(tái)發(fā)動(dòng)機(jī)著陸的算例都是迭代7次收斂。圖8展示了使用單臺(tái)和兩臺(tái)發(fā)動(dòng)機(jī)著陸時(shí)優(yōu)化的控制量，兩個(gè)算例推力變化趨勢(shì)相近，兩臺(tái)發(fā)動(dòng)機(jī)著陸時(shí)大推力持續(xù)時(shí)間更短。單發(fā)動(dòng)機(jī)著陸消耗燃料10.727 t，兩臺(tái)發(fā)動(dòng)機(jī)著陸消耗燃料9.773 t，少消耗954 kg燃料。由此可見(jiàn)飛行器的推力設(shè)計(jì)對(duì)其著陸燃料消耗也是有一定影響的。

圖8 單臺(tái)發(fā)動(dòng)機(jī)和兩臺(tái)發(fā)動(dòng)機(jī)著陸的控制量對(duì)比Fig.8 The differences of controls between one-engine landing and two-engine landing

3.4 先翻轉(zhuǎn)再著陸和邊翻轉(zhuǎn)邊著陸的比較

上文的仿真算例中，目標(biāo)函數(shù)中的懲罰系數(shù)μ都取0.3，當(dāng)目標(biāo)函數(shù)中μ取0時(shí)目標(biāo)函數(shù)就變成了最省燃料，因?yàn)椴豢紤]對(duì)翻轉(zhuǎn)過(guò)程高度差的優(yōu)化，可能出現(xiàn)翻轉(zhuǎn)剛結(jié)束飛行器就接近地面的情況，這里把這種優(yōu)化目標(biāo)對(duì)應(yīng)的優(yōu)化問(wèn)題P0稱(chēng)為邊翻轉(zhuǎn)邊著陸問(wèn)題；μ= 0.3 時(shí)目標(biāo)函數(shù)表示先翻轉(zhuǎn)再著陸的同時(shí)節(jié)省燃料，這種優(yōu)化目標(biāo)對(duì)應(yīng)的問(wèn)題P0就是本文主要討論的先翻轉(zhuǎn)再著陸問(wèn)題。本節(jié)將通過(guò)仿真比較兩種著陸方式的結(jié)果。

先翻轉(zhuǎn)再著陸和邊翻轉(zhuǎn)邊著陸兩個(gè)算例都迭代7 次收斂，邊翻轉(zhuǎn)邊著陸時(shí)燃料消耗8.844 t，比先翻轉(zhuǎn)再著陸少消耗929 kg燃料。圖9展示了兩個(gè)算例的仿真結(jié)果對(duì)比。圖9a是飛行器俯仰角和俯仰角速度變化曲線(xiàn)，可以看出邊翻轉(zhuǎn)邊著陸（μ= 0）時(shí)翻轉(zhuǎn)過(guò)程變長(zhǎng)，翻轉(zhuǎn)過(guò)程變慢。圖9b 顯示了質(zhì)量變化，邊翻轉(zhuǎn)邊著陸更省燃料，因?yàn)槠淠繕?biāo)函數(shù)就是最省燃料。圖9c 和圖9d 是兩個(gè)算例推力和發(fā)動(dòng)機(jī)擺角的對(duì)比。邊翻轉(zhuǎn)邊著陸的前期使用最小推力，因?yàn)椴恍枰焖俜D(zhuǎn)；著陸前先最小推力再最大推力。這個(gè)邊翻轉(zhuǎn)邊著陸的算例中翻轉(zhuǎn)段耗時(shí)較短，有些算例中邊翻轉(zhuǎn)邊著陸問(wèn)題中翻轉(zhuǎn)過(guò)程可以占整個(gè)飛行時(shí)間的一半以上。

圖9 先翻轉(zhuǎn)再著陸和邊翻轉(zhuǎn)邊著陸的結(jié)果對(duì)比Fig.9 The comparison between flipping before landing case and flipping while landing case

4 結(jié) 論

本文對(duì)類(lèi)星艦飛行器翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題進(jìn)行了建模，設(shè)計(jì)了用于求解該問(wèn)題的序列凸優(yōu)化算法。為了提高算法的收斂性，該算法引入虛擬控制，懲罰信賴(lài)域，并給出了較好的初值計(jì)算方法。通過(guò)數(shù)值仿真可以發(fā)現(xiàn)，在求解考慮姿態(tài)動(dòng)力學(xué)的類(lèi)星艦飛行器翻轉(zhuǎn)著陸軌跡優(yōu)化問(wèn)題這樣一個(gè)強(qiáng)非線(xiàn)性軌跡優(yōu)化問(wèn)題時(shí)，與一般的序列凸優(yōu)化算法相比本文設(shè)計(jì)的算法收斂性有很大提高。本文給出的翻轉(zhuǎn)著陸過(guò)程的軌跡優(yōu)化結(jié)果可以為相關(guān)飛行器的研發(fā)提供一定的參考。