基于序列凸優(yōu)化的上面級遠(yuǎn)程交會軌跡優(yōu)化

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

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

0 引言

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

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

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

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

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

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

1.1 問題建模

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

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

圖1 中：軸方向為地心指向上面級初始位置；為上面級繞地心飛過的角度；為上面級位置矢量；為上面級初始軌道半徑；為目標(biāo)器軌道半徑；為推力矢量；為速度矢量；v和v分別為徑向速度分量和切向速度分量。以時間為自變量的傳統(tǒng)動力學(xué)方程為

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

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

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

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

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

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

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

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

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

式中：為發(fā)動機最大推力。

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

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

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

1.2 模型轉(zhuǎn)化

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

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

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

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

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

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

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

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

無損松弛后式（13）改寫為

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

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

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

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

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

2 迭代算法設(shè)計

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

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

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

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

當(dāng)≥1 時，利用對偶內(nèi)點法求解式（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 仿真試驗及結(jié)果分析

本文針對末端時間自由的上面級遠(yuǎn)程交會問題進行仿真，以驗證提出的序列凸優(yōu)化算法的有效性。仿真使用的計算機配置為：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），針對燃料約束下時間最優(yōu)和時間約束下燃料最優(yōu)兩種情況進行仿真。首先令=0，=1，即時間最優(yōu)。仿真結(jié)果如下，最優(yōu)軌跡如圖4 所示。

圖4 時間最優(yōu)交會軌跡（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)移時間為2 426 s。質(zhì)量變化曲線如圖9 所示，可以看到在時間最優(yōu)的情況下，上面級可用燃料被全部耗盡，剩余質(zhì)量1 000 kg。由所求推力序列及上面級初始位置速度進行軌道外推，可得最終位置與目標(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)移時間優(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

下面進行燃料最優(yōu)情況下的仿真，令=1，=0，仿真結(jié)果如下：燃料最優(yōu)、時間最優(yōu)和3 500 s雙脈沖轉(zhuǎn)移三條軌跡的對比圖如圖10 所示，可見燃料最優(yōu)軌跡已經(jīng)接近霍曼轉(zhuǎn)移。推力變化曲線如圖11 所示，總推力仍然為bang-bang 形態(tài)。上面級質(zhì)量變化示意圖如圖12 所示，最終剩余質(zhì)量為2 010 kg。燃料最優(yōu)情況下最終轉(zhuǎn)移時間為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 種情況的對比可以看到，本文所建立模型不局限于時間或燃料的單一優(yōu)化目標(biāo)，可以應(yīng)對多種任務(wù)情況。利用CVX 工具箱求解單次凸優(yōu)化問題耗時約2 s，該計算時間與離散化的節(jié)點數(shù)量有關(guān)，算法需要進行10～20 次迭代求出最終的軌跡。

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

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

4 結(jié)束語

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