郭 杰，相 巖，王 肖，史鵬飛，唐勝景

(1. 北京理工大學(xué)宇航學(xué)院，北京 100081；2. 飛行器動(dòng)力學(xué)與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京 100081；3. 中國(guó)運(yùn)載火箭技術(shù)研究院，北京 100076；4. 北京航天自動(dòng)控制研究所，北京 100854)

0 引 言

近年來(lái)，可重復(fù)使用運(yùn)載器成為航空航天領(lǐng)域的研究熱點(diǎn)。隨著SpaceX公司獵鷹火箭回收和復(fù)用的多次成功，垂直起降重復(fù)使用運(yùn)載火箭受到了廣泛關(guān)注，著陸段制導(dǎo)問(wèn)題是火箭垂直回收的關(guān)鍵性問(wèn)題之一。著陸段約束復(fù)雜、飛行環(huán)境高動(dòng)態(tài)變化，制導(dǎo)算法需要具備在線處理多約束優(yōu)化問(wèn)題的能力。隨著計(jì)算機(jī)性能的提升和優(yōu)化方法的發(fā)展，一些學(xué)者提出了計(jì)算制導(dǎo)的概念，即通過(guò)在線實(shí)時(shí)軌跡優(yōu)化得到制導(dǎo)指令，直接輸入動(dòng)力學(xué)方程實(shí)現(xiàn)閉環(huán)最優(yōu)制導(dǎo)。快速、高精度的在線軌跡優(yōu)化方法是閉環(huán)最優(yōu)制導(dǎo)的基礎(chǔ)。

軌跡優(yōu)化方法一般可分為直接法和間接法。在直接法中，凸優(yōu)化方法求解速度快、具備確定收斂準(zhǔn)則的特點(diǎn)使其在數(shù)值計(jì)算方面具有明顯的優(yōu)勢(shì)。然而大多數(shù)飛行器軌跡優(yōu)化問(wèn)題的動(dòng)力學(xué)和約束都是非凸的，因此應(yīng)用凸優(yōu)化方法的難點(diǎn)在于問(wèn)題的凸化。目前常用的凸化方法主要有序列凸化和無(wú)損凸化兩種。

相較于序列凸化，無(wú)損凸化采用變量替換、約束松弛的方式實(shí)現(xiàn)凸化，無(wú)需參考軌跡和附加的信賴(lài)域約束，具有良好的收斂特性。然而無(wú)損凸化依賴(lài)特殊的問(wèn)題形式，對(duì)于飛行高度較低、大氣稠密的火箭子級(jí)，氣動(dòng)力的影響不可忽略，無(wú)法直接應(yīng)用無(wú)損凸化求解。文獻(xiàn)[13]采用無(wú)損凸化處理推力幅值約束，利用平均速度估計(jì)氣動(dòng)阻力進(jìn)而完成火箭垂直回收著陸段在線制導(dǎo)律設(shè)計(jì)，但是采用平均速度估計(jì)氣動(dòng)力誤差較大，導(dǎo)致軌跡優(yōu)化算法精度和制導(dǎo)算法的魯棒性降低。文獻(xiàn)[14]結(jié)合無(wú)損凸化和同倫方法設(shè)計(jì)了火箭子級(jí)著陸同倫迭代凸化軌跡優(yōu)化方法，但是該算法未考慮推力方向和變化率約束，加重了發(fā)動(dòng)機(jī)系統(tǒng)的工作負(fù)擔(dān)。

現(xiàn)有的無(wú)損凸化方法在確定最優(yōu)終端時(shí)間方面研究較少，多將終端時(shí)間設(shè)置為固定值，無(wú)法保證開(kāi)機(jī)-終端時(shí)刻組合的最優(yōu)性。然而，若直接將作為優(yōu)化變量會(huì)破壞原問(wèn)題的凸性，進(jìn)而影響原問(wèn)題良好的收斂特性。另一方面，對(duì)尋優(yōu)則需要對(duì)包含最優(yōu)解區(qū)間的上下界進(jìn)行多次試探，計(jì)算效率低，求解時(shí)間長(zhǎng)，難以在線應(yīng)用。因此如何在著陸段初始時(shí)刻快速確定最優(yōu)終端時(shí)間成為了在火箭回收問(wèn)題中應(yīng)用無(wú)損凸化方法的主要難點(diǎn)。

偽譜法具有求解精度高、算法成熟等優(yōu)點(diǎn)。因此許多學(xué)者將偽譜離散與凸優(yōu)化相結(jié)合求解軌跡優(yōu)化問(wèn)題，文獻(xiàn)[16]將偽譜離散與無(wú)損凸化結(jié)合研究了火星著陸軌跡優(yōu)化問(wèn)題。文獻(xiàn)[17]將偽譜法與序列凸優(yōu)化結(jié)合，在模型預(yù)測(cè)控制框架下解決了火箭垂直回收制導(dǎo)問(wèn)題。然而上述方法采用全局離散策略，問(wèn)題稀疏性差，求解耗時(shí)較長(zhǎng)，難以在線應(yīng)用。

綜上所述，現(xiàn)有的無(wú)損凸化方法研究中在確定最優(yōu)終端時(shí)間方面研究較少，并且難以實(shí)現(xiàn)軌跡優(yōu)化算法求解精度和計(jì)算速度的有效權(quán)衡。本文面向可重復(fù)使用運(yùn)載火箭子級(jí)回收問(wèn)題，提出了一種帶有最優(yōu)終端時(shí)間估計(jì)策略的hp偽譜同倫凸優(yōu)化在線軌跡規(guī)劃算法。主要貢獻(xiàn)為：1)提出了一種最優(yōu)終端時(shí)間快速估計(jì)策略，能夠快速估計(jì)最優(yōu)終端時(shí)間，同時(shí)保證優(yōu)化結(jié)果的最優(yōu)性；2)所提方法無(wú)需參考軌跡，能夠有效處理氣動(dòng)力和非凸約束，在保證一定求解精度的條件下有效提升了算法的計(jì)算速度。

1 問(wèn)題描述

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

不考慮地球自轉(zhuǎn)，將地球視為均勻圓球，忽略哥氏慣性力的影響，考慮阻力和升力，參考文獻(xiàn)[18]中的發(fā)射坐標(biāo)系定義著陸點(diǎn)坐標(biāo)系如圖1所示。圖中，坐標(biāo)系原點(diǎn)為著陸點(diǎn)，軸垂直地面向上，軸與初始時(shí)刻的速度矢量在地面上投影的方向一致，軸與,軸呈右手坐標(biāo)系。

圖1 著陸點(diǎn)坐標(biāo)系Fig.1 Reference frame of the landing point

在著陸點(diǎn)坐標(biāo)系下建立可重復(fù)使用運(yùn)載火箭子級(jí)著陸段三自由度動(dòng)力學(xué)方程如下

(1)

(2)

式中：和分別為升力系數(shù)和阻力系數(shù)；為火箭子級(jí)參考特征面積；為大氣密度，其計(jì)算公式為

=exp(-)

(3)

式中：為海平面密度；為參考高度；為火箭子級(jí)飛行高度。

對(duì)于火箭垂直回收過(guò)程，令狀態(tài)變量=[,,]，控制變量==[,,]，則系統(tǒng)狀態(tài)空間表達(dá)式為

(4)

1.2 軌跡優(yōu)化問(wèn)題

火箭子級(jí)垂直回收需滿足定點(diǎn)軟著陸的任務(wù)需求，因此定義軌跡優(yōu)化問(wèn)題的初始和終端約束為

(5)

式中：和分別為著陸段初始和終端時(shí)刻；和分別為火箭子級(jí)著陸段初始質(zhì)量和結(jié)構(gòu)干重；對(duì)于運(yùn)載火箭子級(jí)一般有=,=。進(jìn)一步考慮控制量約束，其中推力幅值和方向約束為

(6)

(7)

式中：和為總推力幅值的上下限；max和max為最大推力方向角。

此外，考慮到火箭子級(jí)發(fā)動(dòng)機(jī)推力調(diào)節(jié)能力有限，還需要對(duì)推力變化率施加約束

(8)

另一方面，由于著陸段起始高度較低，火箭子級(jí)飛行速度相對(duì)較小，而空氣密度較大，因此僅考慮動(dòng)壓作為過(guò)程約束

(9)

選擇燃料消耗最少即終端剩余質(zhì)量最大為性能指標(biāo)

=-()

(10)

綜上所述，可重復(fù)使用運(yùn)載火箭子級(jí)著陸段軌跡優(yōu)化問(wèn)題可描述為

問(wèn)題0：

(11)

問(wèn)題0是一個(gè)連續(xù)的終端時(shí)間自由的最優(yōu)控制問(wèn)題，其優(yōu)化變量為=[(),(),]，對(duì)應(yīng)的最優(yōu)控制量為=[(),(),]。

2 帶有最優(yōu)終端時(shí)間估計(jì)的軌跡優(yōu)化算法

2.1 問(wèn)題凸化

無(wú)損凸化

針對(duì)問(wèn)題0中非凸形式的動(dòng)力學(xué)方程和約束，無(wú)法直接應(yīng)用凸優(yōu)化方法求解。本節(jié)主要考慮通過(guò)無(wú)損凸化將非凸問(wèn)題轉(zhuǎn)化為凸問(wèn)題。

首先針對(duì)非凸推力幅值約束，引入變量，則式(6)可以轉(zhuǎn)化為

≤≤

(12)

(13)

基于凸優(yōu)化理論和問(wèn)題的具體形式，對(duì)式(13)進(jìn)行松弛變換，變換后的約束為

(14)

關(guān)于這一凸化變換無(wú)損性的證明可參考文獻(xiàn)[12]和文獻(xiàn)[14]，本文中直接應(yīng)用結(jié)論，則動(dòng)力學(xué)方程可變換為

(15)

基于上述變換，問(wèn)題0中控制變量可增廣變換為=[,]，定義變換后的問(wèn)題為問(wèn)題1，具體形式如下

問(wèn)題1：

(16)

進(jìn)一步對(duì)上述動(dòng)力學(xué)方程中的非線性項(xiàng)進(jìn)行處理，作如下變量替換

(17)

變換后的動(dòng)力學(xué)方程如下

(18)

相應(yīng)地，推力約束和目標(biāo)函數(shù)變換為

·e-≤≤·e-

(19)

(20)

=min-()

(21)

根據(jù)凸優(yōu)化和解析幾何理論，易知式(20)所定義的可行域?yàn)槎A錐，是一種典型的凸約束。

同倫方法

對(duì)于式(18)中的非凸氣動(dòng)力項(xiàng)和式(19)所表示的非凸約束，采用同倫方法(Homotopy approach)對(duì)其進(jìn)行凸化。在數(shù)學(xué)上，兩個(gè)拓?fù)淇臻g如果可以通過(guò)一系列連續(xù)變化從一個(gè)變到另一個(gè)，那么就稱(chēng)這兩個(gè)拓?fù)淇臻g同倫。對(duì)于函數(shù)而言，若存在連續(xù)映射(,)，使得連續(xù)函數(shù)(),()滿足

(22)

則稱(chēng)()與()同倫，(,)為同倫映射，為同倫因子，由上式可知，隨著同倫因子從0變化到1，函數(shù)()逐漸延拓為()。

值得注意的是，應(yīng)用同倫方法的前提條件是同倫過(guò)程中不改變系統(tǒng)的性質(zhì)，下面對(duì)本問(wèn)題中同倫方法的適用性進(jìn)行簡(jiǎn)要分析。

不同于以氣動(dòng)力作為主要控制力的滑翔再入飛行器，垂直回收運(yùn)載火箭的主要控制力為火箭發(fā)動(dòng)機(jī)推力。此外，在垂直回收過(guò)程中火箭子級(jí)飛行速度相對(duì)較小，氣動(dòng)力相對(duì)發(fā)動(dòng)機(jī)推力要小很多，本文采用的仿真參數(shù)對(duì)應(yīng)的作用力曲線如圖2所示。

圖2 飛行過(guò)程中所受作用力幅值Fig.2 Magnitudes of the acting forces during flight

從圖中可以看出，氣動(dòng)力小于推力，且由于推力對(duì)飛行速度的主要控制作用，氣動(dòng)力曲線也在一定程度上受到了推力曲線的影響。

通過(guò)上述分析可知，氣動(dòng)力為火箭子級(jí)的非主要控制力，可以應(yīng)用同倫方法凸化問(wèn)題。同倫凸優(yōu)化算法的具體設(shè)計(jì)過(guò)程如下。

首先不考慮氣動(dòng)力的影響，即在式(18)中刪去氣動(dòng)力加速度項(xiàng)，基于文獻(xiàn)[11]中的處理方式對(duì)式(19)進(jìn)行Taylor展開(kāi)，將其變換為

·e-[(1-(-)+(-)2)]≤

≤·e-[(1-(-)]

(23)

由式(23)可得，不等式左側(cè)為二階錐約束,右側(cè)為線性約束，對(duì)于凸化后的問(wèn)題應(yīng)用原-對(duì)偶內(nèi)點(diǎn)法(Primal-dual interior point method, PDIPM)求解凸優(yōu)化問(wèn)題獲得初始解。隨后，在動(dòng)力學(xué)方程中引入同倫因子，將動(dòng)力學(xué)方程、推力約束以及過(guò)程約束變換為

(24)

·e-≤≤·e-

(25)

(26)

(27)

式中：

需要說(shuō)明的是，上述狀態(tài)空間表達(dá)式中的系數(shù)矩陣與控制矩陣在同倫過(guò)程中為線性時(shí)變剖面形式，具體取值由上一次優(yōu)化結(jié)果決定。

綜上所述，凸化后的問(wèn)題可以描述為

問(wèn)題2：

(28)

2.2 問(wèn)題離散

本文采用flipped Radau偽譜法(flipped Radau pseudospectral method, FRPM)離散連續(xù)最優(yōu)控制問(wèn)題2，將離散后的參數(shù)優(yōu)化問(wèn)題記為問(wèn)題3。FRPM是一種非對(duì)稱(chēng)的多區(qū)間偽譜法，不同于Radau偽譜法(RPM)，F(xiàn)RPM將終端時(shí)間設(shè)置為配點(diǎn)，可以更好地定義終端約束。

在FRPM中，首先定義新的時(shí)間變量，將時(shí)間區(qū)間從[0,]映射到[-1,1]，定義仿射變換關(guān)系如下

(29)

然后選取個(gè)Flipped Legendre Gauss Radau(FLGR)配點(diǎn)以及=-1作為離散點(diǎn)，故在∈[-1,1]上共有個(gè)配點(diǎn)和+1個(gè)離散點(diǎn)。定義階Lagrange多項(xiàng)式()為如下形式

(30)

相應(yīng)地，狀態(tài)量和控制量可以近似為

(31)

對(duì)上式中的狀態(tài)量求導(dǎo)即可得到狀態(tài)量微分的近似形式

(32)

因此動(dòng)力學(xué)方程的離散形式為

=1,2,…,

(33)

式中：∈×(+1)為FRPM微分矩陣。此外，目標(biāo)函數(shù)為Mayer型性能指標(biāo)，其離散形式為

=[(1)]

(34)

相應(yīng)地，初始和終端約束變換為

(35)

推力約束變換為

·e-≤≤·e-,?∈[-1,1]

(36)

(37)

(38)

(39)

過(guò)程約束變換為

(40)

進(jìn)一步地，為加快求解速度、提升求解效率，采用hp偽譜法代替全局偽譜法對(duì)問(wèn)題進(jìn)行離散化處理。本文為簡(jiǎn)化問(wèn)題，僅考慮子區(qū)間個(gè)數(shù)h和配點(diǎn)個(gè)數(shù)p為常值的情況。采用hp-FRPM離散軌跡優(yōu)化問(wèn)題增加了雅各比矩陣的稀疏性，加快了軌跡優(yōu)化問(wèn)題的求解速度。

在本文中用上標(biāo)表示第個(gè)子區(qū)間，下標(biāo)表示子區(qū)間內(nèi)的第個(gè)離散點(diǎn)，則動(dòng)力學(xué)方程式(33)的hp離散形式為

(41)

此外，考慮各子區(qū)間之間的連接條件

(42)

綜上所述，采用hp-FRPM離散后的參數(shù)優(yōu)化問(wèn)題可描述為

問(wèn)題3：

(43)

2.3 最優(yōu)終端時(shí)間快速估計(jì)

現(xiàn)有的無(wú)損凸化方法在確定最優(yōu)終端時(shí)間方面研究較少，多將設(shè)置為固定值，在一定程度上對(duì)求解結(jié)果施加了限制條件，不能完全保證求解結(jié)果的最優(yōu)性。

但是，若直接將納入軌跡優(yōu)化問(wèn)題框架則會(huì)破壞原問(wèn)題良好的收斂特性，難以應(yīng)用無(wú)損凸化求解。另一方面，若采用傳統(tǒng)的線性搜索算法對(duì)尋優(yōu)則最優(yōu)解區(qū)間的上下界難以確定，導(dǎo)致搜索算法收斂速度慢，難以在線應(yīng)用。

基于上述分析，本文基于最優(yōu)控制理論和垂直起降可重復(fù)使用運(yùn)載火箭的運(yùn)動(dòng)學(xué)特性,采用解析方法估計(jì)包含最優(yōu)解區(qū)間的上下界，進(jìn)一步設(shè)計(jì)二次插值型搜索策略完成最優(yōu)終端時(shí)間的快速估計(jì)。

(44)

(45)

式中：為切換系數(shù)。

進(jìn)一步地，由火箭子級(jí)動(dòng)力學(xué)式(1)中第二式可得，不考慮氣動(dòng)力時(shí)火箭子級(jí)加速度幅值為

(46)

雖然火箭子級(jí)加速度幅值為時(shí)變值，難以分析其具體變化趨勢(shì)，但結(jié)合式(45)中的燃料最優(yōu)控制形式和火箭子級(jí)變質(zhì)量飛行特性，可以通過(guò)定性分析式(46)得出，在垂直回收過(guò)程中加速度總體呈增大趨勢(shì)，終端時(shí)刻加速度最大。同時(shí)，根據(jù)終端時(shí)刻火箭子級(jí)姿態(tài)約束和最優(yōu)控制理論可知，終端時(shí)刻火箭子級(jí)推力方向沿著陸點(diǎn)坐標(biāo)系軸方向且幅值為。通過(guò)上述分析可得，火箭子級(jí)終端加速度具體形式為

(47)

對(duì)于上式中的終端質(zhì)量，結(jié)合動(dòng)力學(xué)方程式(1)中的第三式和燃料最優(yōu)控制問(wèn)題解式(45)估計(jì)其下界

(48)

則終端加速度的上界為

(49)

至此，完成終端加速度上界的解析估計(jì)。另一方面，由于火箭子級(jí)初始位置和速度已知，根據(jù)運(yùn)動(dòng)學(xué)關(guān)系，構(gòu)建多項(xiàng)式函數(shù)來(lái)估計(jì)速度，進(jìn)而求解最優(yōu)終端時(shí)間上下界，原理如圖3所示。

圖3 最優(yōu)終端時(shí)間上下界解析多項(xiàng)式估計(jì)Fig.3 Analytic polynomial estimation of upper and lower bounds of optimal terminal time

(50)

(51)

由上式可知，終端時(shí)刻火箭子級(jí)加速度的估計(jì)值為

(52)

進(jìn)一步根據(jù)運(yùn)動(dòng)學(xué)公式，速度曲線對(duì)時(shí)間的積分為位置，即

(53)

將式(50)代入上式求解可得

(54)

(55)

(56)

綜上所述，帶有最優(yōu)終端時(shí)間快速估計(jì)策略的火箭垂直回收hp偽譜同倫凸優(yōu)化算法流程如圖4所示。

圖4 算法流程Fig.4 Flow chart of the algorithm

3 仿真校驗(yàn)

本節(jié)通過(guò)數(shù)值仿真對(duì)提出的帶有最優(yōu)終端時(shí)間快速估計(jì)策略的hp偽譜同倫凸優(yōu)化算法進(jìn)行驗(yàn)證，并將其與現(xiàn)有的軌跡優(yōu)化算法進(jìn)行對(duì)比，分析該軌跡優(yōu)化算法的求解精度、計(jì)算時(shí)間以及收斂性能。

本文使用的仿真參數(shù)及過(guò)程約束見(jiàn)表1和表2，火箭子級(jí)初始位置和速度設(shè)置為=[-2700,9800,0]m、=[180,-400,0]m/s。所有程序的仿真環(huán)境均為Intel Core i7-7700HQ、2.80 GHz和Windows 10操作系統(tǒng)的筆記本電腦，仿真軟件為Matlab，對(duì)于生成的凸問(wèn)題調(diào)用Mosek 9.0工具包采用原始對(duì)偶內(nèi)點(diǎn)法進(jìn)行求解。

3.1 最優(yōu)終端時(shí)間快速估計(jì)策略仿真驗(yàn)證

表1 火箭子級(jí)仿真參數(shù)Table 1 Simulation parameters of the launch vehicle stages

表2 約束邊界值Table 2 Boundary values of constraints

由于偽譜法將終端時(shí)間納入優(yōu)化框架，通過(guò)求解軌跡優(yōu)化問(wèn)題來(lái)獲得最優(yōu)終端時(shí)間，因此可以將偽譜法的求解結(jié)果作為基準(zhǔn)。從表3中可以看出，四種方法的求解結(jié)果相差不大。進(jìn)一步對(duì)比求解耗時(shí)，偽譜法求解時(shí)間明顯較長(zhǎng)。通過(guò)對(duì)比文獻(xiàn)[21]方法和改進(jìn)文獻(xiàn)[21]方法的求解時(shí)間可以得出，采用本文的上下界估計(jì)方法可以有效縮小估計(jì)區(qū)間，加快計(jì)算速度，提升計(jì)算效率。本文方法的求解時(shí)間為1.58 s，遠(yuǎn)小于其余三種方法的求解時(shí)間，說(shuō)明了采用二次插值法具有超線性收斂的特性，相較于基于線性收斂準(zhǔn)則的二分法，收斂速度大幅提升。同時(shí)本文方法的求解結(jié)果相較于基準(zhǔn)結(jié)果誤差較小，實(shí)現(xiàn)了具備一定精度的快速估計(jì)，適合在線應(yīng)用。

表3 算法性能對(duì)比Table 3 Performance comparison of algorithms

3.2 算法仿真分析及最優(yōu)性驗(yàn)證

本節(jié)仿真以國(guó)內(nèi)外廣泛應(yīng)用的精度高、適用性強(qiáng)的GPOPS Ⅱ軟件包的求解結(jié)果作為基準(zhǔn)，對(duì)提出的hp偽譜同倫凸優(yōu)化算法即hp-PHCP算法的最優(yōu)性進(jìn)行驗(yàn)證。基于GPOPS Ⅱ軟件包的求解程序采用hp網(wǎng)格自適應(yīng)更新策略，調(diào)用經(jīng)典的稀疏非線性求解器(SNOPT)求解轉(zhuǎn)化后的非線性規(guī)劃問(wèn)題。hp-PHCP算法中總離散點(diǎn)個(gè)數(shù)為81，其中子區(qū)間個(gè)數(shù)=4，配點(diǎn)個(gè)數(shù)=20，飛行終端時(shí)間采用3.1節(jié)中的優(yōu)化結(jié)果，同倫參數(shù)設(shè)置為=04。上述兩個(gè)算法的迭代終止條件參數(shù)均為=10。為表示方便，在圖例中將利用GPOPS Ⅱ軟件包求解記為原算法，將hp-PHCP算法記為本文算法，仿真結(jié)果如圖5至圖13所示。

圖5 hp-PHCP求解SOCP問(wèn)題時(shí)間Fig.5 Solving time of SOCP problem using hp-PHCP

圖6 動(dòng)壓-時(shí)間曲線Fig.6 Curve of dynamic pressure-time

圖7 推力擺角-時(shí)間曲線Fig.7 Curve of thrust swing angle-time

圖8 推力幅值迭代過(guò)程Fig.8 Iteration history of thrust amplitude

圖9 位置-時(shí)間曲線Fig.9 Curve of position-time

圖10 速度-時(shí)間曲線Fig.10 Curve of velocity-time

圖11 質(zhì)量-時(shí)間曲線Fig.11 Curve of mass-time

圖12 推力分量-時(shí)間曲線Fig.12 Curve of thrust component-time

圖13 推力幅值-時(shí)間曲線Fig.13 Curve of thrust amplitude-time

根據(jù)仿真結(jié)果，hp-PHCP算法共迭代6次完成收斂，求解二階錐規(guī)劃(SOCP)問(wèn)題總時(shí)間為1.011 s，每次求解時(shí)間如圖5所示。從圖中可以看出，hp-PHCP算法求解SOCP問(wèn)題的平均時(shí)間約為0.169 s，最大求解時(shí)間為0.347 s，出現(xiàn)在求解迭代啟動(dòng)解中，這是由于啟動(dòng)解求解過(guò)程無(wú)需參考軌跡，采用全局尋優(yōu)造成了時(shí)間較長(zhǎng)。從圖6至圖8可以看出，求解結(jié)果滿足動(dòng)壓約束、推力幅值和方向約束。特別地，從圖8的局部放大圖可知，在前3次迭代中推力切換時(shí)間跳變較大，且推力幅值約束并不滿足。這是由于進(jìn)行了控制量和質(zhì)量的變量替換，而氣動(dòng)力尚未完全加入動(dòng)力學(xué)系統(tǒng)，在同倫過(guò)程中以上次優(yōu)化結(jié)果為基準(zhǔn)，但質(zhì)量剖面與氣動(dòng)力剖面的近似精度不足，導(dǎo)致變換后的推力幅值約束近似精度下降。這一問(wèn)題在氣動(dòng)力完全引入動(dòng)力學(xué)方程中得到了顯著改善，從圖8的局部放大圖可以看出，在第四次迭代后，同倫參數(shù)=1即完成氣動(dòng)力引入，總推力曲線迅速收斂，推力幅值約束得到了滿足。從圖9至圖13可以看出，hp-PHCP算法與GPOPS優(yōu)化軟件包求得的狀態(tài)變量曲線基本重合，控制變量曲線相差不大，因此可以說(shuō)明，提出的hp-PHCP算法具備一定的局部最優(yōu)性。特別地，從圖12可以看出，由于推力方向和變化率約束的存在，分推力曲線并不平滑，然而狀態(tài)變量曲線平滑，這表明了hp-PHCP算法在嚴(yán)格滿足控制約束的同時(shí)實(shí)現(xiàn)了平穩(wěn)飛行。此外，從圖13還可以看出，松弛變量與總推力曲線完全重合，進(jìn)一步驗(yàn)證了凸化過(guò)程的無(wú)損性。

3.3 算法性能對(duì)比分析

為進(jìn)一步分析本文提出hp-PHCP算法的性能，采用對(duì)比分析的方法對(duì)三種不同離散方法的軌跡優(yōu)化算法進(jìn)行數(shù)值仿真。第一種是采用基于梯形法離散的同倫凸優(yōu)化算法(THCP)，第二種是基于全局flipped Radau偽譜法離散的同倫凸優(yōu)化算法(PHCP)，第三種是本文提出的基于hp flipped Radau偽譜法離散的同倫凸優(yōu)化算法(hp-PHCP)。三種算法的總離散點(diǎn)數(shù)均為81，飛行終端時(shí)間采用3.1節(jié)中的優(yōu)化結(jié)果，同倫參數(shù)設(shè)置為=04。其中hp-PHCP算法中子區(qū)間個(gè)數(shù)=2，配點(diǎn)個(gè)數(shù)=40。仿真結(jié)果如表4和圖14至圖18所示。

表4 三種算法優(yōu)化結(jié)果對(duì)比Table 4 Optimization results comparison of three algorithms

圖14 位置-時(shí)間曲線Fig.14 Curve of position-time

圖15 速度-時(shí)間曲線Fig.15 Curve of velocity-time

圖16 質(zhì)量-時(shí)間曲線Fig.16 Curve of mass-time

圖17 推力分量-時(shí)間曲線Fig.17 Curve of thrust component-time

圖18 推力幅值-時(shí)間曲線Fig.18 Curve of thrust amplitude-time

從圖14至圖18可以看出，三種軌跡優(yōu)化算法求解得到的狀態(tài)變量和控制變量曲線基本重合，總推力曲線呈現(xiàn)出明顯的Bang-Bang特性，且切換過(guò)程平滑無(wú)吉布斯(Gibbs)現(xiàn)象。此外，三種軌跡優(yōu)化算法均在第6次迭代時(shí)完成收斂，收斂性能良好，收斂速度快。

表4給出了三種軌跡優(yōu)化算法優(yōu)化結(jié)果對(duì)比，作為以在線應(yīng)用為目標(biāo)的軌跡優(yōu)化算法，求解精度和計(jì)算速度是衡量算法優(yōu)劣的重要標(biāo)準(zhǔn)。三種算法中THCP求解時(shí)間最短，但高度和速度的終端誤差大；而PHCP算法雖然計(jì)算精度最高，但其求解時(shí)間長(zhǎng)。相較于THCP和PHCP，本文提出的hp-PHCP算法終端誤差與PHCP算法為同一數(shù)量級(jí)，但是求解時(shí)間不到PHCP算法的一半，在保證一定求解精度的條件下有效提升了算法的計(jì)算速度，使其具備了在線應(yīng)用的潛力。

4 結(jié) 論

面向可重復(fù)使用運(yùn)載火箭子級(jí)回收問(wèn)題，提出了一種帶有最優(yōu)終端時(shí)間估計(jì)策略的hp偽譜同倫凸優(yōu)化在線軌跡規(guī)劃算法。取得的主要成果如下：

1) 最優(yōu)終端時(shí)間快速估計(jì)策略能夠快速估計(jì)最優(yōu)終端時(shí)間，保證了開(kāi)機(jī)-終端時(shí)刻組合的最優(yōu)性，同時(shí)能夠最大限度實(shí)現(xiàn)燃料消耗最小，進(jìn)一步提升運(yùn)載火箭的經(jīng)濟(jì)效益。

2) hp偽譜同倫凸優(yōu)化軌跡優(yōu)化算法(hp-PHCP)能夠有效處理非凸約束，無(wú)需參考軌跡，分推力指令平滑連續(xù)，收斂性能良好，收斂速度快，能夠保證求解結(jié)果的最優(yōu)性。

3) 相較于THCP算法和PHCP算法，hp-PHCP算法基于hp多階段偽譜離散策略，進(jìn)一步提升了軌跡優(yōu)化問(wèn)題的稀疏性，在保證求解精度的同時(shí)有效縮短了計(jì)算時(shí)間，具備在線應(yīng)用潛力。