王訪寒，周如好，余薛浩，黃 飛，陳海朋

（上海航天控制技術(shù)研究所，上海 201109）

0 引言

上面級(jí)是航天運(yùn)輸系統(tǒng)的重要組成部分之一。與基礎(chǔ)級(jí)火箭相比，上面級(jí)具有多次啟動(dòng)、長(zhǎng)期在軌、自主飛行、軌道機(jī)動(dòng)能力強(qiáng)等特點(diǎn)，工作時(shí)通常已經(jīng)進(jìn)入地球軌道，具備較強(qiáng)的靈活性和通用性，可以適應(yīng)多種不同的任務(wù)情況。近年來(lái)，隨著在軌服務(wù)需求的不斷增加，上面級(jí)將繼續(xù)向著自主化、智能化等方向發(fā)展，其中自主交會(huì)和軌跡優(yōu)化是其中的關(guān)鍵技術(shù)，也是研究的熱點(diǎn)問題。

目前常見的軌道優(yōu)化方法有采用龐特里亞金極小值原理的間接法、以高斯偽譜法為代表的直接法、以遺傳算法為代表的智能算法等。間接法將軌跡優(yōu)化問題轉(zhuǎn)化為哈密頓邊值問題，但是存在推導(dǎo)過程復(fù)雜、協(xié)態(tài)變量初始值難以估計(jì)等問題。在此基礎(chǔ)上發(fā)展出的直接法將問題轉(zhuǎn)化為非線性規(guī)劃問題，對(duì)初始猜測(cè)不敏感，但是計(jì)算量較大，常用于離線軌跡規(guī)劃。還有智能算法同樣具有計(jì)算耗時(shí)較長(zhǎng)、實(shí)時(shí)性差等問題。

在諸多軌道優(yōu)化方法中，凸優(yōu)化方法在實(shí)時(shí)性方面具有非常明顯的優(yōu)勢(shì)。尤其是對(duì)偶內(nèi)點(diǎn)法的發(fā)現(xiàn)，可以確保二階錐規(guī)劃問題在有限步迭代內(nèi)收斂到最優(yōu)解。ACIKMESE 等將凸優(yōu)化技術(shù)成功運(yùn)用于火星登陸器軟著陸的最優(yōu)控制問題。LIU 和LU 等采用無(wú)損凸化的方法將交會(huì)對(duì)接中的非凸問題轉(zhuǎn)化為凸問題并求出優(yōu)化軌跡。WANG 等采用角度而非時(shí)間作為自變量，求解了小推力航天器軌道轉(zhuǎn)移問題。林曉輝等基于凸優(yōu)化研究了月球精確定點(diǎn)軟著陸問題。李鑫等基于序列凸優(yōu)化求解了固定時(shí)間遠(yuǎn)程交會(huì)問題。池賢彬等運(yùn)用凸優(yōu)化制導(dǎo)技術(shù)研究了近程自主交會(huì)問題。

目前，相關(guān)文獻(xiàn)中還有一些待解決的問題。首先，研究大多是基于固定時(shí)間假設(shè)求解燃料最優(yōu)問題，并未說(shuō)明如何確定軌道轉(zhuǎn)移所用的時(shí)間，在沒有時(shí)間限制的常規(guī)交會(huì)任務(wù)中，如果給定的時(shí)間過長(zhǎng)或過短，則優(yōu)化出的軌跡并不是真正意義的燃料最優(yōu)；其次，在空間救援或目標(biāo)監(jiān)視等緊急的特殊任務(wù)中，會(huì)要求上面級(jí)在燃料約束下用最短時(shí)間與目標(biāo)器完成交會(huì)，此時(shí)固定時(shí)間的軌跡優(yōu)化方法不再適用；最后，現(xiàn)有文獻(xiàn)對(duì)迭代初始值的猜測(cè)缺乏詳細(xì)說(shuō)明，一般是采用初位置與末位置的線性插值作為位置的初始猜測(cè)，雖然序列凸優(yōu)化具有對(duì)初值不敏感的優(yōu)點(diǎn)，但是在上面級(jí)大范圍軌道轉(zhuǎn)移的情況下，這種偏差較大的初值在某些情況下會(huì)對(duì)收斂性產(chǎn)生一定影響，所以需要一個(gè)較為接近實(shí)際情況的初值來(lái)保證收斂性，加快收斂速度。

針對(duì)以上問題，本文采用地心距代替時(shí)間作為自變量，建立了末端時(shí)間自由的上面級(jí)遠(yuǎn)程交會(huì)模型，并采用Lambert 雙脈沖變軌的計(jì)算結(jié)果作為迭代初值，最后在時(shí)間最優(yōu)和燃料最優(yōu)兩種任務(wù)情況下進(jìn)行了仿真驗(yàn)證，求出了兩種情況下的最優(yōu)軌跡和推力控制序列。

1 基于凸優(yōu)化的上面級(jí)遠(yuǎn)程交會(huì)模型

1.1 問題建模

本文針對(duì)末端時(shí)間自由的共面圓軌道遠(yuǎn)程交會(huì)問題進(jìn)行研究，為便于問題描述，建立了極坐標(biāo)系，如圖1 所示。

圖1 軌道轉(zhuǎn)移極坐標(biāo)系示意圖Fig.1 Polar coordinates for orbital transfer

圖1 中：軸方向?yàn)榈匦闹赶蛏厦婕?jí)初始位置；為上面級(jí)繞地心飛過的角度；為上面級(jí)位置矢量；為上面級(jí)初始軌道半徑；為目標(biāo)器軌道半徑；為推力矢量；為速度矢量；v和v分別為徑向速度分量和切向速度分量。以時(shí)間為自變量的傳統(tǒng)動(dòng)力學(xué)方程為

式中：為上面級(jí)到地心的距離；為發(fā)動(dòng)機(jī)總推力大??；F和F分別為發(fā)動(dòng)機(jī)總推力在徑向和切向的分量；為地球引力常數(shù)；為發(fā)動(dòng)機(jī)噴氣速度。

以時(shí)間作為自變量常用于固定時(shí)間下燃料最優(yōu)問題，而對(duì)于末端時(shí)間自由的遠(yuǎn)程交會(huì)問題，一種做法是內(nèi)層采用固定時(shí)間的傳統(tǒng)優(yōu)化方法，外層采用二分法或智能算法對(duì)最優(yōu)時(shí)間進(jìn)行搜索，這種方法沒有真正意義上將時(shí)間作為優(yōu)化變量加入到模型中，每一次外層搜索都進(jìn)行了一次固定時(shí)間的軌跡優(yōu)化，會(huì)大大增加總計(jì)算時(shí)間，無(wú)法應(yīng)用于線上計(jì)算。

另一種做法是利用歸一化將時(shí)間變換到區(qū)間[0，1]，再當(dāng)成固定時(shí)間問題進(jìn)行求解。然而這種方法會(huì)導(dǎo)致末端時(shí)間這一未知量始終存在，后續(xù)的問題凸化將變得十分困難。同理，轉(zhuǎn)移角度在末端時(shí)間自由的交會(huì)問題中也不適合作為自變量。造成這一問題的根源為：當(dāng)末端時(shí)間不固定時(shí)，時(shí)間和轉(zhuǎn)移角度都沒有一個(gè)確定的上界，模型離散化后待求向量的維數(shù)無(wú)法確定，導(dǎo)致無(wú)法求解。

為解決上述問題，本文采用上面級(jí)到地心的距離作為自變量，在徑向速度不小于零的前提下，的取值范圍為[，]且單調(diào)，把時(shí)間作為待優(yōu)化變量加入到模型中，利用時(shí)間與地心距的關(guān)系改寫動(dòng)力學(xué)方程如下：

應(yīng)滿足的初始條件為

式中：為上面級(jí)初始質(zhì)量。

應(yīng)滿足的末端條件為

式中：為目標(biāo)器初始相位角；為遠(yuǎn)程交會(huì)可用燃料質(zhì)量。此外，還可以根據(jù)任務(wù)需求增加時(shí)間約束。

基于有限推力的控制約束為

式中：為發(fā)動(dòng)機(jī)最大推力。

優(yōu)化指標(biāo)可表示為

式中：和分別為交會(huì)時(shí)間和燃料消耗的加權(quán)因子；Δ和Δ分別為軌道轉(zhuǎn)移所用時(shí)間和消耗燃料質(zhì)量。當(dāng)=0 且=1 時(shí)，燃料最優(yōu)；當(dāng)=1 且=0 時(shí)，時(shí)間最優(yōu)。

此外，還應(yīng)根據(jù)實(shí)際任務(wù)需要，將變量約束在一個(gè)大致合理的范圍內(nèi)，同時(shí)用小量代替約束條件中徑向速度為0 的情況，以避免出現(xiàn)奇異問題。

1.2 模型轉(zhuǎn)化

若想應(yīng)用凸優(yōu)化解決該問題，就需要將原始問題中的非凸部分進(jìn)行凸化。上述模型中的非凸項(xiàng)主要存在于動(dòng)力學(xué)中的速度位置、質(zhì)量變化和推力約束。

首先，為消除質(zhì)量變化產(chǎn)生的非凸部分，定義新變量如下：

對(duì)式（9）兩邊同時(shí)對(duì)半徑求導(dǎo)得

由式（10）可以看到，用變量代替質(zhì)量，并用代替了/項(xiàng)，消除了原模型中質(zhì)量變化造成的非線性部分。

其次，為消除推力約束產(chǎn)生的非凸部分，定義新變量：

則式（5）和式（6）所表示有限推力約束可寫為

式（13）～式（14）中表示的推力約束在三維空間中是一個(gè)二階錐的表面，是非凸的，通過無(wú)損松弛將其轉(zhuǎn)化成凸約束，在三維空間中表示實(shí)心二階錐，如圖2 所示。

圖2 控制約束無(wú)損凸化Fig.2 Lossless convexification of control constraint

無(wú)損松弛后式（13）改寫為

將式（14）右側(cè)在參考軌跡z處泰勒展開，取一階近似得

最后，通過序列線性化的方法來(lái)處理動(dòng)力學(xué)中剩余的非凸部分，經(jīng)過1.1 節(jié)中的變量替換后，動(dòng)力學(xué)方程的形式為

式 中：=[，，v，v，]為狀態(tài)變量；=[u，u，]為控制變量。向量()和矩陣()的表達(dá)式為

假設(shè)上一次迭代后產(chǎn)生軌跡的狀態(tài)變量為，則動(dòng)力學(xué)方程線性化后可近似為

至此已將原模型中所有非凸部分凸化，記第次迭代的狀態(tài)向量為，使用一階差分進(jìn)行離散化處理，步長(zhǎng)即采樣距離設(shè)為R，該問題轉(zhuǎn)化后的最終形式為

2 迭代算法設(shè)計(jì)

進(jìn)行序列凸優(yōu)化的迭代求解需要初始參考軌跡，一般是采用初末位置線性插值的方法給出初始軌跡。雖然序列凸優(yōu)化具有對(duì)初值不敏感的優(yōu)點(diǎn)，但是一個(gè)合理的初值可以保證收斂性，減少迭代次數(shù)，增加算法實(shí)時(shí)性。本文結(jié)合上面級(jí)發(fā)動(dòng)機(jī)推力較大的特點(diǎn)，采用Lambert 雙脈沖轉(zhuǎn)移軌道作為初始參考。由于Lambert 法已有成熟的求解步驟，在固定時(shí)間條件下可快速求解出轉(zhuǎn)移軌道，相關(guān)文獻(xiàn)已有詳細(xì)研究，此處不再贅述。

求解序列凸優(yōu)化的算法步驟如下：

給定脈沖軌道轉(zhuǎn)移時(shí)間，利用Lambert法計(jì)算出脈沖轉(zhuǎn)移軌道作為初始猜測(cè)。

初始化：令=0，根據(jù)初始猜測(cè)的轉(zhuǎn)移軌道，計(jì)算出，進(jìn)而求出()、()、()的值。

當(dāng)≥1 時(shí)，利用對(duì)偶內(nèi)點(diǎn)法求解式（23）～式（24）表示的凸優(yōu)化問題，得到優(yōu)化后的軌跡和控制向量。

判斷所求軌跡是否滿足收斂條件，收斂條件如下：

式中：為允許誤差。若滿足條件則執(zhí)行步驟6，若不滿足條件則執(zhí)行步驟5。

利用求得的軌跡更新()、()、()的值，并執(zhí)行步驟3。

所求的最優(yōu)軌跡和控制向量為和。

算法流程如圖3 所示。

圖3 序列凸優(yōu)化算法流程Fig.3 Flow chart of the sequential convex optimization algorithm

3 仿真試驗(yàn)及結(jié)果分析

本文針對(duì)末端時(shí)間自由的上面級(jí)遠(yuǎn)程交會(huì)問題進(jìn)行仿真，以驗(yàn)證提出的序列凸優(yōu)化算法的有效性。仿真使用的計(jì)算機(jī)配置為：CPU i5-7500@3.4 GHz，4 GB 內(nèi)存，在Matlab 環(huán)境中使用CVX 工具箱求解迭代中的單次凸優(yōu)化問題。

表1 仿真參數(shù)表Tab.1 Simulation parameters

優(yōu)化目標(biāo)參考式（23），針對(duì)燃料約束下時(shí)間最優(yōu)和時(shí)間約束下燃料最優(yōu)兩種情況進(jìn)行仿真。首先令=0，=1，即時(shí)間最優(yōu)。仿真結(jié)果如下，最優(yōu)軌跡如圖4 所示。

圖4 時(shí)間最優(yōu)交會(huì)軌跡（k1=0，k2=1）Fig.4 Minimum time rendezvous trajectory when k1=0 and k2=1

推力變化曲線如圖5～圖7 所示，可以看出總推力呈現(xiàn)典型的bang-bang 形態(tài)。目標(biāo)函數(shù)的優(yōu)化過程如圖8 所示，算法經(jīng)過12 次迭代后收斂，最終轉(zhuǎn)移時(shí)間為2 426 s。質(zhì)量變化曲線如圖9 所示，可以看到在時(shí)間最優(yōu)的情況下，上面級(jí)可用燃料被全部耗盡，剩余質(zhì)量1 000 kg。由所求推力序列及上面級(jí)初始位置速度進(jìn)行軌道外推，可得最終位置與目標(biāo)器相差149 km，為最終軌道半徑的1.23%，證明了該算法的有效性。

圖5 總推力變化（k1=0，k2=1）Fig.5 Variation of the total thrust when k1=0 and k2=1

圖6 切向推力變化（k1=0，k2=1）Fig.6 Variation of the tangential thrust when k1=0 and k2=1

圖7 徑向推力變化（k1=0，k2=1）Fig.7 Variation of the radial thrust when k1=0 and k2=1

圖8 軌道轉(zhuǎn)移時(shí)間優(yōu)化過程（k1=0，k2=1）Fig.8 Optimization process of the orbital transfer time when k1=0 and k2=1

圖9 質(zhì)量變化（k1=0，k2=1）Fig.9 Variation of the mass when k1=0 and k2=1

下面進(jìn)行燃料最優(yōu)情況下的仿真，令=1，=0，仿真結(jié)果如下：燃料最優(yōu)、時(shí)間最優(yōu)和3 500 s雙脈沖轉(zhuǎn)移三條軌跡的對(duì)比圖如圖10 所示，可見燃料最優(yōu)軌跡已經(jīng)接近霍曼轉(zhuǎn)移。推力變化曲線如圖11 所示，總推力仍然為bang-bang 形態(tài)。上面級(jí)質(zhì)量變化示意圖如圖12 所示，最終剩余質(zhì)量為2 010 kg。燃料最優(yōu)情況下最終轉(zhuǎn)移時(shí)間為4 735 s。

圖10 3 種優(yōu)化軌跡示意圖（k1=1，k2=0）Fig.10 Three optimized trajectories when k1=1 and k2=0

圖11 總推力變化（k1=1，k2=0）Fig.11 Variation of the total thrust when k1=1 and k2=0

圖12 質(zhì)量變化（k1=1，k2=0）Fig.12 Variation of the mass when k1=1 and k2=0

通過2 種情況的對(duì)比可以看到，本文所建立模型不局限于時(shí)間或燃料的單一優(yōu)化目標(biāo)，可以應(yīng)對(duì)多種任務(wù)情況。利用CVX 工具箱求解單次凸優(yōu)化問題耗時(shí)約2 s，該計(jì)算時(shí)間與離散化的節(jié)點(diǎn)數(shù)量有關(guān)，算法需要進(jìn)行10～20 次迭代求出最終的軌跡。

將本文提出的序列凸優(yōu)化算法與遺傳算法進(jìn)行對(duì)比，見表2。序列凸優(yōu)化算法的耗時(shí)為28 s，遺傳算法的優(yōu)化速度為305 s。文獻(xiàn)［20-21］表明偽譜法在軌跡優(yōu)化中的計(jì)算耗時(shí)一般為分鐘級(jí)，本文提出的序列凸優(yōu)化相比于偽譜法和遺傳算法具有明顯的速度優(yōu)勢(shì)，具備線上計(jì)算的潛力。

表2 仿真結(jié)果對(duì)比Tab.2 Comparison of the simulation results

4 結(jié)束語(yǔ)

本文研究了基于序列凸優(yōu)化的上面級(jí)遠(yuǎn)程交會(huì)軌跡優(yōu)化問題。針對(duì)末端時(shí)間自由的遠(yuǎn)程交會(huì)問題，采用地心距代替時(shí)間作為自變量，通過無(wú)損凸化將問題轉(zhuǎn)化為凸問題，使用Lambert 脈沖轉(zhuǎn)移軌道作為迭代初值，分別求出了時(shí)間最優(yōu)和燃料最優(yōu)任務(wù)情況下的轉(zhuǎn)移軌跡和推力序列。如果針對(duì)軌道轉(zhuǎn)移問題開發(fā)定制化的求解器，單次求解凸優(yōu)化問題的時(shí)間還能進(jìn)一步縮小，未來(lái)可應(yīng)用于在線計(jì)算。但本文未考慮攝動(dòng)因素，并且模型局限于圓軌道之間的交會(huì)任務(wù)，后續(xù)研究可考慮將模型擴(kuò)展至異面橢圓軌道的交會(huì)任務(wù)中。