基于混合優(yōu)化算法的RLV 覆蓋區(qū)求解

李惠峰 孫國慶 何睿智

（北京航空航天大學宇航學院，北京 100191）

1 引言

近年來由于可重復使用運載器（Reusable Launch Vehicle，RLV）擁有大的縱程和橫程機動能力，受到各國科研人員廣泛的關注，各個國家均已投入大量人力、物力進行研究。當飛行器再入時，飛行器需要較長一段時間在大氣層中飛行，將受到熱流密度、動壓、過載的嚴格約束，又由于飛行器需要進行主動姿態(tài)控制，給飛行器材料和內部設施熱絕緣帶來巨大的挑戰(zhàn)，所以飛行器安全、可靠的再入成為非常重要的問題。

再入飛行器覆蓋區(qū)是指飛行器以某一初始條件再入在地球表面能夠獲得著陸的范圍，是評估再入飛行器航程覆蓋能力的重要性能指標。當再入飛行器返回時，通過覆蓋區(qū)的信息和飛行器的狀態(tài)來確定飛行器的著陸點。

國內外的研究工作中主要有四種覆蓋區(qū)求解方法：通過坐標旋轉將覆蓋區(qū)求解轉化為在自由縱程下最大橫程的參數(shù)搜索問題［1］；利用直接法轉換為軌跡優(yōu)化問題［2－3］；通過最大和最小阻力包線近似判定覆蓋區(qū)外邊界和內邊界［4］；文獻［5－6］中論述在規(guī)定的縱程下，最近軌跡和最大橫程問題，兩個問題解決方法相同，將覆蓋區(qū)問題轉化為最近軌跡下的單一參數(shù)搜索問題。大部分方法在求解覆蓋區(qū)問題時都需要進行參數(shù)搜索，而且往往是多參數(shù)搜索，傳統(tǒng)的參數(shù)優(yōu)化主要有基于極大值原理的間接法和基于非線性規(guī)劃理論的直接法［7］。但無論是直接法還是間接法，都需要根據(jù)工程經驗設置參數(shù)初值，然后通過試錯法對這些參數(shù)進行調整，而目標函數(shù)對待優(yōu)化參數(shù)的微小變化十分敏感，給全局尋優(yōu)帶來很大困難，很難找到全局最優(yōu)解，并要花費大量的時間。所以優(yōu)化結果的好壞在很大程度上取決于某些優(yōu)化變量的初始猜測值以及優(yōu)化算法的性能。

本文針對覆蓋區(qū)求解過程中優(yōu)化參數(shù)的初值選取和優(yōu)化算法性能與效率問題，利用多學科優(yōu)化軟件，找到參數(shù)優(yōu)化的初值點，基于參數(shù)的響應面特性，在分析全局優(yōu)化算法和局部優(yōu)化算法優(yōu)缺點的基礎上，設計了混合優(yōu)化算法，將具有較好全局優(yōu)化性能的遺傳算法和局部優(yōu)化性能突出的模式搜索法相結合，仿真結果表明混合優(yōu)化算法具有精度高、收斂快等特點，能夠準確且高效地獲得再入覆蓋區(qū)，有很好的工程實用價值。

圖1 覆蓋區(qū)示意Fig.1 Illustration of footprint

2 覆蓋區(qū)問題建模

（1）覆蓋區(qū)問題描述

本文在求解覆蓋區(qū)時，采用定縱程下的最大橫程方法來求解。如圖1所示，E為再入點，ED 是最大的縱程，AD 和CD 是給定的縱程（EB1，EB2，…，EBK）下最大橫程的集合，AC 是最小縱程的集合，覆蓋區(qū)定義為被AD、CD 和AC 連接包圍的范圍。這里不考慮最小縱程的集合AC，主要考慮外邊界AD 和CD 的求解。

（2）無動力再入運動方程

考慮地球為旋轉圓球旋轉時，引起的科里奧利力（Coriolis Force）和牽引力引起的影響，RLV 無動力再入飛行器，在極坐標系下的量綱一運動方程為

為了提高數(shù)值計算精度，對各變量進行歸一化處理，即用原變量除以對應的歸一化系數(shù)，其中r為地心矢徑長度，歸一化系數(shù)為地球半徑R0；V 是速度，歸一化系數(shù)為L，D 是氣動升力和阻力，歸一化系數(shù)為mgn，相應的升力系數(shù)為CL，阻力系數(shù)為CD，m 是飛行器質量；γ是航跡傾角；σ是滾轉角；θ，φ是經度和緯度；ψ 是方位角；Ω 是地球自轉角速度。時間t采用進行歸一化，控制變量為滾轉角σ和攻角α。如無特別說明，下文中變量均為歸一化變量。

采用以下假設條件：重力加速度為常值gn；滿足準平衡滑翔條件；高度為常值且r＝1；不考慮地球自轉，忽略哥氏加速度與牽連加速度。經簡化處理得到三自由度運動方程：

（3）平衡滑翔約束條件

對于再入飛行器而言，再入過程受到動力學和物理學約束，比如氣動加熱，過載及動壓約束［6］，再入飛行器以高超聲速在大氣層中飛行時，空氣受到強烈的壓縮和摩擦作用，大部分能量轉化成熱能，導致周圍空氣溫度急劇升高，而且熱能迅速向飛行器表面?zhèn)鬟f，引起嚴重的氣動加熱問題。如果氣動加熱過度，會損壞飛行器和內部設備，所以對溫度的限制十分有必要；而機體結構和運載器表面放熱材料的強度決定了飛行器的動壓約束；飛行器在飛行時受到載荷作用，當載荷超過承受極限，飛行器結構和內部儀表可能會遭到破壞，因此也需要對飛行器所受的總過載進行限制。這些常規(guī)約束對再入過程影響很大，此外為保證平滑再入飛行，還要考慮準平衡滑翔約束。

對具有高升阻比的再入飛行器，在進行滑翔飛行時，其法向加速度較小，從而導致航跡傾角以及航跡傾角變化率很小。由此，高聲速飛行器三自由度再入運動方程組中，表示航跡傾角變化的微分方程，就可以通過將航跡傾角變化率近似為零進行簡化，這個代數(shù)式就是平衡滑翔條件。平衡滑翔條件給出了飛行器滑翔軌跡上速度與高度很重要的關系式，通過將其應用到升力體飛行器再入相關問題的研究中，能夠顯著降低問題的維數(shù)與難度。但是，在對平衡滑翔的實際應用中發(fā)現(xiàn)，飛行器的航跡傾角變化率不會理想地保持為零，它隨時間在很小的范圍內緩慢變化，很多情況下航跡傾角出現(xiàn)小幅值的長周期震蕩動態(tài)。在這種情況下，平衡滑翔條件在實際應用中改稱為準平衡滑翔條件（Quasi－Equilibrium Glide Condition，QEGC）：

利用QEGC中給出的速度－高度關系可以將飛行器再入中的過程約束轉化為對滾轉角控制量σ最大值的約束，熱流率約束、動壓約束以及過載約束對應的滾轉角最大值分別用σmax＿Q（V）、σmax＿q（V）和σmax＿n（V）表示，具體值計算：

式中 Qmax、qmax和nmax為飛行器所允許的最大熱流率值、最大動壓值以及最大過載值；KQ為與飛行器自身有關的常數(shù)；KL＝R0SrefCL／2m，Sref為飛行器參考面積，m為飛行器質量。綜上所述，滾轉角的最大取值為

3 最優(yōu)控制問題

（1）給定縱程下的最大橫程問題

對定縱程下最大橫程問題，通過軌跡最優(yōu)控制，使得飛行器從再入點E（由經緯度θE和φE表示），在給定縱程SP（對應點P，由經緯度θP和φP表示）下，達到最大橫程SF（對應點F，由經緯度θF和φF表示）。優(yōu)化目標是使得橫程SF最大，性能泛函

終點約束條件主要有兩個：達到預設的末端能量；EP 與PF 的夾角為直角：

式中 rf和Vf是終端高度和速度。

哈密頓函數(shù)可寫為

由于時間t并不顯含于哈密頓函數(shù)中，所以有：

式中 P為伴隨變量矢量；x為狀態(tài)變量矢量；c0為常數(shù)。則可得：

同時由H2和H3的表達式可以得到：

經過一系列推導（具體過程參見文獻 ［8］），最終可得到滾轉角最優(yōu)控制律為

在滾轉角表達式中，分子分母都是積分常數(shù)c1，c2以及c3的線性表達式，由此可見，滾轉角表達式的值只與其中兩個獨立變量有關，最終的優(yōu)化目標也就是要找到這兩個獨立變量的具體值。結合最大滾轉角約束條件可得：

式中 σ＊為設計的滾轉角控制律，且｜σ｜＜π／2。由于末端航跡偏角不固定，可以是任意值，由最優(yōu)化理論易知：

參數(shù)c1，c2和c3中只有兩個獨立變量，假設：c2＝sink，c3＝cosk，上式整理可得：

至此定縱程下最大橫程優(yōu)化問題表述為：在滿足終點約束條件的前提下，利用一系列的尋優(yōu)算法求解最優(yōu)控制參數(shù)c1和k的值。

（2）最大縱程問題

對于最大縱程問題，性能泛函為

終點約束條件為末端能量約束，如式（1）所示，對于第2節(jié)中的再入動力學方程，從初始再入點開始，采用零度滾轉角控制律對軌跡進行積分，當再入軌跡的能量變量滿足終端能量停機條件時，軌跡積分停止，得到的軌跡終點即為最大縱程點D，計算得到該軌跡對應的縱程即為最大縱程ED。

4 參數(shù)優(yōu)化分析

在第3節(jié)的問題求解中把能量條件設置為停機條件，當飛行器能量因子達到預設值，那一時刻的飛行器所處點就作為再入軌跡飛行的終端點，需要求解兩個參數(shù)（c1，k）以使其滿足非線性終端約束方程：

對于該雙參數(shù)尋優(yōu)問題，基于這兩個非線性函數(shù)，構造一個性能泛函J（c1，k）：

令ω1＝ω2＝1，得到性能泛函J（c1，k）在區(qū)間c1∈［－1.5，1.5］，k∈［－π，π］的掃描圖（圖2）。這是一個很不規(guī)則的參數(shù)響應面，具有很多局部凹槽，并且對參數(shù)的微小變化十分敏感?？梢钥闯鰧⒎蔷€性函數(shù)J1（c1，k）與J2（c1，k）取平方和會造成新的性能泛函J（c1，k）出現(xiàn)很多凹凸面，使得在基于梯度方法進行搜索時很容易陷入局部最優(yōu)的情況，很難找到全局最優(yōu)解。

圖2 性能泛函J 掃描Fig.2 Performance index J

通過對整個可行域內的參數(shù)響應分析，可以發(fā)現(xiàn)對于該具體的飛行器尋優(yōu)問題，在整個可行域內有4個最優(yōu)點，通過編程求解方程必定可以得到初始的一個解，對于其他的解，則可以通過一定的波動或者改變搜索范圍進行查找。但由于滾轉角控制律中，分子分母均為參數(shù)的線性表達式，所以導致了這4個最優(yōu)解實際的物理意義有所不同，實際上只有兩個解是最優(yōu)解，一個是左邊最優(yōu)點，一個是右邊最優(yōu)點。另外兩個“最優(yōu)點”雖然對于該飛行器再入階段的滾轉角控制量是一樣的，但在二階特性中，是不滿足條件的，應該舍棄。

優(yōu)化算法可以分為兩大類：全局優(yōu)化算法、局部優(yōu)化算法。相比較而言，全局優(yōu)化算法收斂速度較慢，而局部優(yōu)化算法收斂速度快。但局部優(yōu)化算法存在一個初始點的選取問題，如果初始點選取不夠準確，將有可能使搜索進入局部極值，結果不準確。

基于多學科優(yōu)化軟件OPTIMUS所集成的優(yōu)化算法對典型的優(yōu)化算法進行測試。全局優(yōu)化算法包括：差分進化、自適應進化、模擬退火。局部優(yōu)化算法包括：序列二次規(guī)劃、非線性二次規(guī)劃、廣義簡約梯度。六種優(yōu)化算法所對應的優(yōu)化結果精度分別為：0.000 26，0.000 22，0.000 38，2.556 13×10－5，0.000 56，3.391 79×10－5。通過對比分析優(yōu)化結果，得出如下結論：1）全局優(yōu)化算法不需要給定初始點，在整個參數(shù)可行域內進行概率搜索，得到的優(yōu)化結果往往較接近于全局最優(yōu)點，但精度不高；2）如果給出的初始點位置較好，局部算法要明顯短于全局優(yōu)化算法，而且局部算法可以設置非常短的迭代步長，精度較高；3）局部算法出錯的可能性較高，這是由局部算法本身的特性所導致的，局部搜索算法依賴于給定的初始點，敏感性非常強，如果給定的初始點離全局最優(yōu)點較遠的話，往往容易陷入局部極值點，得到錯誤的結果。

從平臺全局算法的優(yōu)化結果，可以看出當?shù)玫阶顑?yōu)結果的同時，也得到了全局最優(yōu)點，可以將其作為初值點寫入局部優(yōu)化算法中，以縮短局部尋優(yōu)時間。

5 混合優(yōu)化方法

根據(jù)第4節(jié)分析，提出一種結合不同方式的優(yōu)缺點的混合優(yōu)化方法，將全局優(yōu)化得到的優(yōu)化結果點作為局部優(yōu)化的初始點進行進一步精確搜索，得到更加精確解，滿足全局最優(yōu)性，而且可以極大地降低錯誤率。

5.1 遺傳算法

遺傳算法的特點在于其直接對結構對象進行操作，不需要導數(shù)信息或者要求函數(shù)連續(xù)，具有很好的全局搜索能力，不需要給出初始點，采用概率化的尋優(yōu)方式能自適應地調整搜索方向。其算法分為以下步驟：

1）第一步：編碼。將可行域中的優(yōu)化設計向量進行編碼。在規(guī)定的搜索區(qū)間內（c1可行域［－1，1］，k可行域 ［－π，π］），對兩個變量進行編碼，編碼形式為碼長100的二進制碼，前50個表示參數(shù)c1，后50個表示參數(shù)k。

2）第二步：形成初始群體。隨機選擇若干個體形成初始群體，初始種群個數(shù)設置為1 000個。

3）第三步：計算適應值。與目標函數(shù)直接相關，有相同的極值點。

4）第四步：選擇。將適應值累加，利用賭輪盤法進行選擇。

5）第五步：交叉。選擇交叉概率為0.6，從交配池中隨機取出幾個個體配對，確定交叉點，將選中的個體在交叉點以后的編碼數(shù)字進行交換。

6）第六步：變異。選定變異概率為0.01，對應編碼的每一位也用隨機小數(shù)表示，小于變異概率就求補運算，進行變異。

對于初始種群中的所有樣本點，計算其適應值，利用賭輪盤的方法進行選擇交叉，同時以之前設定好的變異概率對種群中的“基因”進行變異。不斷迭代循環(huán)，每一步記錄下該步迭代的最優(yōu)點，直至最終的遺傳代數(shù)達到設定值20。搜索結束，記錄下迭代最后的“最優(yōu)點”。

5.2 模式搜索法

模式搜索法是一種從幾何意義上提出的參數(shù)直接尋優(yōu)局部搜索算法。該算法不需要對目標函數(shù)求導，可以解決不可導函數(shù)以及求導復雜函數(shù)的參數(shù)尋優(yōu)?；痉椒椋簭某跏键c開始，依次向四周進行搜索，一旦發(fā)現(xiàn)比當前點的函數(shù)值更優(yōu)的點，就將該點作為下一個搜索周期的起始點進行下一步迭代，直至最后的點的四周均沒有比當前點更優(yōu)的點時，尋優(yōu)結束。模式搜索法受初始搜索點的選取影響很大，具有較快的收斂速度，但并不具有全局最優(yōu)特性，容易進入局部最優(yōu)“凹槽”，搜索步長也需要根據(jù)具體問題進行選取。具體步驟為：

1）確定初始搜索點。將遺傳算法得到的“最優(yōu)點”作為模式搜索法局部搜索的初始點。

2）求周圍點函數(shù)值。以初始搜索點為中心，建立一個立方體，將該立方體的八個頂點作為第一次迭代的樣本點，分別計算其函數(shù)值。

3）判斷是否搜索到更優(yōu)點。將上一步得到的八個樣本點函數(shù)值進行統(tǒng)計，如果其中有一點的函數(shù)值比中心初始搜索點的函數(shù)值更優(yōu)，則記為找到下一個搜索初始點，重新回到步驟2）進行下一步搜索。否則，進行下一步。

4）改變搜索步長。增大搜索步長，在更大的搜索范圍內進行搜索，如果找到更優(yōu)點，回到步驟2）。如果還沒有找到更優(yōu)點，進一步增大搜索步長，直至搜索步長超過一定值，停止搜索。

通過這樣一種全局與局部優(yōu)化算法組合的方法，結合了兩種優(yōu)化方式的共同優(yōu)點，兼顧了全局性和精確性，其精度能達到10－5。

5.3 覆蓋區(qū)求解策略

在解決給定縱程下的最大橫程問題后，為了求出完整的覆蓋區(qū)，需要首先找到最大縱程，然后在最大縱程軌跡上選取不同的縱程點，分別應用給定縱程下的最大橫程求解方法，結合混合優(yōu)化算法進行優(yōu)化搜索，從而求得最大橫程以及所對應的終端點，將所有的終端點連起來，就獲得了再入覆蓋區(qū)。由于混合優(yōu)化算法在快速性和準確性上的優(yōu)勢，即使選取較多的縱程點，也不會影響覆蓋區(qū)的求解效率，從而可以快速地獲得飛行器的精確覆蓋區(qū)。

6 仿真結果

飛行器基本參數(shù)：質量37 362.9kg，參考面積149.388m2，最大允許熱流率為794 425W／m2，最大允許動壓為14 364Pa，最大允許過載為6 gn，攻角α變化范圍為 ［0°，45°］；馬赫數(shù)Ma 變化范圍為 ［3，25］；高度變化范圍為 ［0km，120km］。最大升阻比在1～1.4之間，隨馬赫數(shù)的增加，升阻比減??；在正攻角情況下，升阻比隨攻角增大先增大后減??；最大升阻比出現(xiàn)在8°～10°攻角附近。仿真初始條件和終端條件如表1所示。

表1 仿真初始條件與終端條件Tab.1 Initial and terminal conditions of simulation

先給出在給定縱程為9 871.1km 時的最大橫程問題的仿真結果，給定縱程下的對應左右兩邊各有一個最大橫程，且由于受地球自轉的影響，右邊橫程大于左邊橫程，圖3是左右兩條軌跡的最優(yōu)滾轉角控制曲線。由于存在滾轉角約束，在超過最大約束值時，取最大約束值；在速度5 000～8 000m／s左右，約束起作用，如圖4所示。圖5是再入軌跡，很好地滿足了熱流、動壓和過載的約束，圖6是航跡傾角。

圖3 滾轉角控制曲線Fig.3 Bank angle profile

圖4 存在過程約束時的滾轉角曲線Fig.4 Bank angle profile with path constraint

圖5 速度－高度軌跡曲線Fig.5 Velocity－altitude trajectory

圖6 航跡傾角變化曲線Fig.6 Flight path angle profile

圖7 覆蓋區(qū)Fig.7 Footprint

在最大縱程上取多個點（這里取6 個點，如果要進一步提高覆蓋區(qū)的精度，可以選取更多的點），然后對每一個縱程求最大橫程，把每個橫程的終點連起來，就得到了飛行器的覆蓋區(qū)，如圖7 所示?？梢钥吹?，右半部分要比左半部分大，這是由地球自轉所導致的。

［1］VINH N X.Optimal trajectories in atmospheric flight［M］.New York：Elsevier Scientific Publishing Company，1981.

［2］AMITABH SARAF，JAMES A LEAVITT，KENNETH D MEASE.Landing footprint computation for entry vehicles［C］.AIAA Guidance，Navigation，and Control Conference and Exhibit，Providence，Rhode Island，2004.

［3］雍恩米，唐國金，陳磊.高超聲速無動力遠程滑翔飛行器多約束條件下的軌跡快速生成 ［J］.宇航學報，2008，29（1）：46－52.YONG ENMI，TANG GUOJIN，CHEN LEI.Rapid trajectory planning for hypersonic unpowered long－range reentry vehicles with multi－constraints［J］.Journal of Astronautics，2008，29（1）：46－52.

［4］HALE JR N W，LAMOTTE N O，GARNER T W.Operational experience with hypersonic flight of the space shuttle：Proceedings of 2002 AIAA Guidance，Navigation，and Control（GN＆C）Conference，AIAA，Orlando，F(xiàn)L，2002 ［C］.

［5］SHEN Z J.On－board three－dimensional constrained entry flight trajectory generation ［D］.Ames：Iowa State University，2002.

［6］LU P，XUE S B.Rapid generation of accurate entry landing footprints［J］.Journal of Guidance，Control and Dynamics，2010，33（3）：756－767.

［7］T B JOHN.Survey of numerical methods for trajectory optimization ［J］.Journal of Guidance，Control and Dynamics，1998，21（3）：193－207.

［8］LU P.Entry trajectory optimization with analytical feedback bank angle law ［C］.AIAA Guidance，Navigation and Control Conference and Exhibit，Honolulu，2008.