基于序列凸優化的上面級遠程交會軌跡優化

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

（上海航天控制技術研究所，上海 201109）

0 引言

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

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

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

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

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

1 基于凸優化的上面級遠程交會模型

1.1 問題建模

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

圖1 軌道轉移極坐標系示意圖Fig.1 Polar coordinates for orbital transfer

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

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

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

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

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

應滿足的初始條件為

式中：為上面級初始質量。

應滿足的末端條件為

式中：為目標器初始相位角；為遠程交會可用燃料質量。此外，還可以根據任務需求增加時間約束。

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

式中：為發動機最大推力。

優化指標可表示為

式中：和分別為交會時間和燃料消耗的加權因子；Δ和Δ分別為軌道轉移所用時間和消耗燃料質量。當=0 且=1 時，燃料最優；當=1 且=0 時，時間最優。

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

1.2 模型轉化

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

首先，為消除質量變化產生的非凸部分，定義新變量如下：

對式（9）兩邊同時對半徑求導得

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

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

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

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

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

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

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

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

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

假設上一次迭代后產生軌跡的狀態變量為，則動力學方程線性化后可近似為

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

2 迭代算法設計

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

求解序列凸優化的算法步驟如下：

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

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

當≥1 時，利用對偶內點法求解式（23）～式（24）表示的凸優化問題，得到優化后的軌跡和控制向量。

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

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

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

所求的最優軌跡和控制向量為和。

算法流程如圖3 所示。

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

3 仿真試驗及結果分析

本文針對末端時間自由的上面級遠程交會問題進行仿真，以驗證提出的序列凸優化算法的有效性。仿真使用的計算機配置為：CPU i5-7500@3.4 GHz，4 GB 內存，在Matlab 環境中使用CVX 工具箱求解迭代中的單次凸優化問題。

表1 仿真參數表Tab.1 Simulation parameters

優化目標參考式（23），針對燃料約束下時間最優和時間約束下燃料最優兩種情況進行仿真。首先令=0，=1，即時間最優。仿真結果如下，最優軌跡如圖4 所示。

圖4 時間最優交會軌跡（k1=0，k2=1）Fig.4 Minimum time rendezvous trajectory when k1=0 and k2=1

推力變化曲線如圖5～圖7 所示，可以看出總推力呈現典型的bang-bang 形態。目標函數的優化過程如圖8 所示，算法經過12 次迭代后收斂，最終轉移時間為2 426 s。質量變化曲線如圖9 所示，可以看到在時間最優的情況下，上面級可用燃料被全部耗盡，剩余質量1 000 kg。由所求推力序列及上面級初始位置速度進行軌道外推，可得最終位置與目標器相差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 軌道轉移時間優化過程（k1=0，k2=1）Fig.8 Optimization process of the orbital transfer time when k1=0 and k2=1

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

下面進行燃料最優情況下的仿真，令=1，=0，仿真結果如下：燃料最優、時間最優和3 500 s雙脈沖轉移三條軌跡的對比圖如圖10 所示，可見燃料最優軌跡已經接近霍曼轉移。推力變化曲線如圖11 所示，總推力仍然為bang-bang 形態。上面級質量變化示意圖如圖12 所示，最終剩余質量為2 010 kg。燃料最優情況下最終轉移時間為4 735 s。

圖10 3 種優化軌跡示意圖（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 質量變化（k1=1，k2=0）Fig.12 Variation of the mass when k1=1 and k2=0

通過2 種情況的對比可以看到，本文所建立模型不局限于時間或燃料的單一優化目標，可以應對多種任務情況。利用CVX 工具箱求解單次凸優化問題耗時約2 s，該計算時間與離散化的節點數量有關，算法需要進行10～20 次迭代求出最終的軌跡。

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

表2 仿真結果對比Tab.2 Comparison of the simulation results

4 結束語

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